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ABSTRACT 


This report deals with progress made on the Grant NSG-3048 during the 
twelve month period beginning October 1, 1981 and ending September 30, 1982. 
The NASA Technical Officer for this period was Dr. Kurt Seldner of Lewis 
Research Center. The director of the research at the University of Notre 
Dame was Dr. Michael K. Sain, who has been assisted by Mr. Stephen Yurkovich, 
a fellow in the Department of Electrical Engineering, by Mr. Joe P. Hill, 
and by Mr. Thomas A. Klingler, research assistants, in the Department of 
Electrical Engineering. Mr. Hill received the degree of Master of Science 
during this period, for his June 1982 thesis entitled "Solution of Non- 
linear Optimal Cor.trol Problems Using the Algebraic Tensor: An Example". 

Mr. Klingler expects to complete rasearch investigations for the Master 
of Science degree very shortly. Mr. Yurkovich may complete requirements 
for the degree of Doctor of Philosophy in 1983. 

Researches during the preceding calendar year have centered on basic 
topics in the modeling and feedback control of nonlinear dynamical systems. 

Of special interest have been the following topics: (1) the development 

of models of tensor type for a digital simulation of the QCSE gas turbine 
engine; (2) the extension, to nonlinear multivariable control system de- 
sign, of the cc-.cepts of total synthesis which trace their roots back to 
certain early investigations under this grant; (3) the role of series de- 
scriptions as they relate to questions of scheduling in the control of 
gas curbine engines; (4) the development of computer-aided design soft- 
x^are for tensor modeling calculations; (5) further enhancement of the 
softwares for linear total synthesis, mentioned above; and (6) calcula- 
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tion cf the first known examples using tensors for nonlinear feedback con- 
trol. 

A number of major milestones have occurred during this year of study. 
Most crucial has been the steady progress of the computer software needed 
to perform tensor modeling and simulation. The advance of this code is 
now making it possible to begin a more systematic examination of tensor 
model identification and order reduction. Increasing availability of 
this capability has made it easier to determine certain of the key trade- 
offs involved with use of tensor models * as for example their increased 
dynamical quality versus their useful region. In what may be the most 
significant theoretical development of the year’s activity, work is under- 
way to evaluate the effects of redefining the groups upon which, and into 
which, nonlinear maps act. For interesting cases, the groups can be re- 
defined in such a manner that the nonlinear maps become linear. Among 
the results already following from this discovery are the definition of 
a variety of new nonlinear sensitivity functions, of the comparison type. 
The same idea, in a different application, has permitted the definition 
of a nonlinear feedback synthesis problem. Finally, a complete calcula- 
tion has been carried out for the feedback tensors in a nonlinear control 
example. This calculation has been of the utmost importance in setting 
goals for the type of software which will be needed for general feedback 
control, with tensors. Preliminary steps to plan such software have 
been set in motion. 
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I. INTRODUCTION 


In this report, we discuss progress which has been made on NASA Grant 
NSG-3048, entitled "Alternatives for Jet Engine Control", during the twelve 
month period beginning on October 1, 1981 and ending on September 30, 1982. 

1.1. Overview of Report 

The report is organized broadly into a main body, consisting of six 
sections, and seven appendices. The principal portion of the main body 
involves Sections II, III, and IV, which deal with the explicit tensor 
calculations involved in a substantial example of nonlinear feedback con- 
trol, based upon tensor system models developed in an earlier year of work 
under this grant. These sections are supported by Section 1.2, which 
gives certain literature background, by Section VI, which contains the 
references, and by Appendix G, which illustrates, in an introductory way, 
typical calculations which become involved. 

Insofar as we know, this substantial example is the first of its 
type to be completed. More comment on the insights which follow from 
the main body of the report are provided in Section V, where additional 
discussion is directed to the appendices. 

Appendix E, "A Computer-Aided Design Package for Nonlinear Model 
Applications", is a report on the crucial software developments which 
have been, and concii-ue to be, the backbone of studies in tensor modeling 

CP* 

and simulation. Steady progress in this regard is now making possible 
further work in tensor model identification and order reduction. A com- 
plete description of this software is planned for the next grant period. 
Appendix B, "An Application of Tensor Ideas to Nonlinear Modeling of a 
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Turbofan Jet Engine", illustrates the potential of this developing software 
for work with engine simulations, such as the QCSEE. 

Appendix D, "Controller Schedul^ng: A Possible Algebraic Viewpoint", 

continues the investigation of algebi' ic frameworks which may capture the 
essence of the practical work in control schedules. It is believed that 
these studies may hold one of the keys to reducing gaps between the theory 
of nonlinear control systems and its application. In particular, it is 
planned to use ideas which grow out of this work as a guide to resolving 
the tradeoffs between increased dynamical quality in nonlinear models and 
their region of validity. More study is needed in this area. 

Se tions C and F deal with material on the Total Synthesis Problem 
of multivariable control. In the linear case, this problem traces a part 
of its early roots back to studies supported by this grant. It is a 
problem of feedback synthesis, which has now developed quite a biblio- 
graphy. Appendix F gives a list of that bibliography, as of the date of 
the First American Control Conference last year. Item 18 in Appendix F 
refers to the software associated with this effort, which is ongoing. 

Though not yet ready for distribution, the software has resulted in 
several requests for copies, and it is hoped that limited distribution 
might not be an event too far in the future. Appendix C, "Nonlinear 
Multivariable Design by Total Synthesis", is a part of the effort to 
extend to the nonlinear case. Of particular interest in Appendix C is 
the modification of group structure and scalar multiplication structure 
on the vector spaces of inputs and outputs for certain nonlinear systems. 

It is an amazing fact that, in interesting, nontrivial case, such modi- 
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fications can result in the system becoming linear. The idea is akin to 
the choice of a special coordinate system which fits the geometry of a 
physical problem. We have this idea under intense study, and it seems 
to hold promise in problems <•: order reduction and identification. If 
the work proceeds as planned, more results in this regard should become 
available in next year's report. 

Appendix A contains a chronological bibliography of work undertaken 
in relation to the grant. 

With these introductory remarks, we turn next to the main body of 
the report. 
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1.2 Remarks on Multilinear Feedback 


One subject of this report is the application of concepts from tensor 
algebra to the generation of optimal feedback controllers for nonlinear 
dynamical systems. The primary motivation for a study of this topic came 
from the results presented in [1], where thr local theoretical problem 
was essentially solved. The Importance of this work should be emphasized. 
Although there were works that were previously existln. *,hat used multi- 
linear algebra to study series solutions for a response of a nonlinear 
system, this was the first application of the ideas in tensor algebra to 
an optimal control problem. Thus, we find that the topics of multilinear 
algebra and optimization had been extensively studied, but independently 
of one another. 

The research project has received numerous benefits from this par- 
ticular study. First, we have completed a rather tho "jugh example study 
of the methods in [1]. We believe that this example study is the only 
one of its kind in the literature to date. The details of the example are 
given in Section IV. Second, we have carried out the study without the 
explicit use of the methods of dual spaces and symmetric tensor algebra. 

The purpose of doing so is to gain insight concerning the exact role of 
these two ideas in the nonlinear feedback problem. This was a revealing 
experience; and much insight has been gained. Third, we have an initial 
step in constructing software for such feedback calculations. The im- 
portance of such softwares can scarcely be overestimated. 

The methods of optimal control have, in certain cases, interesting 
relationships with the methods of stability theory. Because of this, we 
have reason to believe that these results may assist in model region design. 
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The treatment of optimal control problems in the history of the lit- 
erature has been quite extensive, with much emphasis placed on the so 
called Linear Quadratic or LQ problem. Broadly speaking, an LQ control 
formulation consists of a finite-dimensional linear discrete - or con- 
tinuous-time dynamic system which is to be controlled in such a way as to 
minimize the value of a performance criterion which is the integral, or 
sum, of quadratic functions of the system state and control variables 
plus a quadratic function of the state at some terminal time, This 

concept provides the following system description: 


x(t) *= A(t) x(t) + B(t) u(t) , 

J = M(x(t^)) f ^ L(x,u,t) dt , 


T T 

L(x,u,t) = X Q(t)x + u R(t)u . 
M(x(t^)) = x^(t^)M x(t^) . 


( 1 . 1 ) 

( 1 . 2 ) 


(1.3a) 

(1.3b) 


Here, we assume that M is symmetric and positive definite, R(t) is 
symmetric and positive definite, and Q(t) is S 3 rmmetric and positive 
semldefinite. So, given the linear system in (1.1) and the cost functional 
in (1.2) satisfying the symmetry and definiteness requirements, we wish 
to find the optimal control, that is, the control which will drive the 
system so as to minimize the cost functional. We do not go into detail 
concerning conditions for existence and uniqueness. Basically, the solu- 
tion of the state regulation problem leads to an optimal feedback system 
with the property that the components of the state vector x(t) are 
kept near zero without excessive expenditure of control energy, which is, 
in essence, the minimization of the cost functional. Thus, given an 
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initial perturbation which is usually in the form of the initial state 
vector x(0), we find that an optimal control should drive the state 
vectors to zero while simultaneously minimizing the selected perforaavice 
index. The notion of desiring the state vector near zero arises from 
the fact that the state variable is defined as an error term [1] which 
measures deviation from the global trajectory. The fact that this per- 
turbation is required to be sufficiently small allows for the traditional 
Taylor's series expansion form. It is well known from the literature 
that we may construct this optimal control as 

u(t) « K^(t) x(t) , (1.4) 

where 

K^(t) = - R"^(t) B^(t) V(t) . (1.5) 

'';.ie matrix V(t) is the solution to the well known matrix Riccati Equa- 
tion 

V(t) = -V(t) A(t) - A^(t) V(t) + V(t) B(t) R"^(t) B'“'^(t) V(t) - Q(t) (1.6) 
which can be readily solved on the digital computer by integrating nu- 
merically backwards in time from the boundary condition 

V(t^) = M . (1.7) 


The problem that we wish to consider in this report is more complex, 
although we will be able to treat the LQ problem as a special case. 
Namely, we wish to assume that we are given a nonlinear system rather 
than a linear system and examine the methods necessary to generate higher 
order controller terras. This control problem is commonly described in 
the literature as being approximately optimal, or sub-optimal. When 
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considering the synthesis of a control function for a typical nonlinear 
system, we will find that an infinite number of terms would be necessary 
in order to formulate a true optimal controller, that is, we must trun- 
cate to a finite number of terms and obtain a performance index that is 
non-minimal. The amount of work that has been done in the generation of 
suboptimal control algorithms is quite extensive. Therefore, in order 
to present an account of the research that has been done, we must be 
succinct and mention only a few of these works. One of the first attempts 
to formulate a description of the optimal control problem as applied to 
nonlinear systems is found in [2], where extensive use is made of the 
methods of Lyapunov and Chetaev. Here, a formal recursive procedure is 
developed to construct a suboptimal control as a function of a power 
series in the states. The work of Lukes [3] extended the concepts intro- 
duced by Al’brekht and provided useful results pertaining to the exis- 
tence and uniqueness of an optimal feedback controller. In [A], we find 
useful applications of the Hamilton-Jacobi-Bellraan approach to a number 
of illustrative examples, as well as a comparison of various techniques 
that can be used to generate higher order controllers. Additional re- 
sults and examples of this method were provided in [5] with particular 
emphasis on the convergence of the procedure. Also, a method for esti- 
mating the degradation in performance caused by the truncation of terms 
in the controller series was presented. We also note the results pre- 
sented in [6], where linear, second order, and third order controller 
expressions were produced using the methods presented in [7]. The im- 
proved response obtained whenever higher order controller terms are 
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considered provides sufficient, motivation to seek another method by which 
these terms may be obtained. This method, of course, is the tensor alge- 
bra. A few of the ideas surrounding the general topic of algebraic system 
theory have previously been applied to problems in optimality, system mod- 
eling, and multivariable feedback loop closures. The latter point was ex- 
tensively studied in [8], where emphasis was placed on the exterior or skew- 
symmetric algebra. The most pertinent work was of course [1] which, as was 
previc'.oly mentioned, provided much of the motivation for this report. An- 
other application of the ideas of tensor algebra to systems problems was in 
[9], which relied heavily on the series expansion concept expressed via the 
tensor product and applied to the basic system description. Particular em- 
phasis was placed on the subject of nonlinear system modeling, with examp- 
les of both homogeneous and nonhomogeneous modeling. These methods proved 
to be quite effective, yielding much improvement over the standard linear 
approximations that are typically used in a nonlinear system for modeling 
purposes. So basically, there would seem to be strong motivation for the 
use of modern algebra in nonlinear systems and control problems. Much of 
what we actually do as systems engineers evolves from modern algebraic con- 
cepts [10]. Since the typical engineer has had little or no exposure to 
the concepts of modern algebra, we will outline these concepts as they are 
related to the solution of our problem. As was noted in [8], modern alge- 
bra frequently provides sufficient algebraic framework within which to ob- 
tain solutions to systems problems with considerably less effort as com- 
pared to other conventional methods. We will find that the tensor alge- 
bra vleT^qDoint provides a useful means of expanding a nonlinear system in 


tenna of vectors and .nntricea, Vv^lrlch are expressed via the tensor product. 


In Section II we provide some useful mathematical background, with 
particular emphasis on the topic of tensor algebra. We examine the pro- 
perties of multilinear mappings and also provide a brief introduction to 
the sym me tr ie tensor algebra structure, whicli will be useful when dif- 
ferentiating the tensor product. In Section IIT we present certain sys- 
tems concepts, which are bosically in the form of series expansions, and 
then derive the necessary results for the generation of the optimal con- 
trol terms. The problem formulation for the LQ optimization can be rec- 
ognized as a special case of the equations wlu're higher order terms are 
included. Finally, Section IV provides tlie application of the results 
derived iu the previous chapter to a formulated example problem. 
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II. SKLKCTKl) ALGKBRAIG BACKGROUND 


The basic purpose of this section Is to provide the reader with the 
necessary concepts from the subject of tensor algebra, which is the main 
vehicle that is used to analyze the control problem. It is realized that 
most readers have not had previous dealings with the somewhat theoretical 
concept of tensor algebra; therefore the treatment of this subject will 
not assume any previous knowledge of the topic. The main feature of the 
algebraic tensor involves the way that it gives ground on dimensionality 
in order to gain the powerful advantage of linearity. 

We begin this section the basic concepts of multilinear mappings 
and the properties that allow us to express these multilinear mappings 
in terms of linear mappings and tensor products. Next, we examine the 
tensor product of linear mappings and develop the associated Kronecker 
product and a few of its properties, which will be useful in later der- 
ivations. The next section deals with the symmetric tensor algebra, 
which will be of great importance whenever we consider the possibilities 
of differentiating the tensor pfiducC, a topic that is considered in the 
final portion of this section. Specifically, wo examine the partial 
derivative problem as related to the tensor product of various copies 
of the state variable x and the control function u. This problem is 
inherently related to the minimization operation that will be studied 
in Section III. 

2.1 Multilinear Mappin.:;s [1,10,11] 

Since many of the ideas surrounding tensor algebra are based on the 
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theory of bilinear and multilinear mappings, we begin with a get\eral def- 
inition of these mappings. Generally speaking, a multilinear function is 
a function of vectors that is linear with respect to each vector variable 
when the others are held constant. This means that 

^ , . . . , 6 y £»•••» (x^ , . . . , , . . . , X^) "b 

3 ij/(Xj^, . . . • • • » * ( 2 . 1 ) 

where and U are vector spaces over a field R; x^, y^ e V^; a , 6 R; 

and i|^:V,xV^x...xv ->-U. The multilinear mapping ip has an image 

1 Z m 

which is not, in general, a subspace of U. A simple counterexample is 
offered in order to illustrate this point. 


We let = V2 be two dimensional spaces with the basis {e^,e2) 
and U be a four dimensional space vjith the basis {f^,f2,f2>f^}* Let 
x^e^^ + X262 and y = ^]_ ^2' 

define a bilinear mapping by 

♦ (x.y) - + V2*4 ’ 


X 


Any vector in U 


u = t c f 
i=l ^ 


is in the image of if and only if it satisfies 

rank 

If we pick u^ = 2 f^ + 2f2 + and U2 = f^ + f^, then obviously 

both of these vectors are in the image of ’i.'. However, subtracting these 
two vectors gives 

“1 ■ “2 “ '1 + ^*2 + ^ 



for which 
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rank 


1 2 
0 1 


= 1 


so that the image is not closed under addition. Since the closure law 
has been violated, we can conclude that Im ip does not form a. linear 
subspace of U. The smallest subspace of U which contains the image 
of ip is called the subspace generated by the image of ip and is de- 
noted by <Im ip> . This subspace is shown by the dotted contour in Figure 
2.1, which depicts the situation for the bilinear case. This minimal sub- 
space becomes a space of tensors when the bilinear map tp is a tensor 
product . 


A particular subset (m=2) of the set of multilinear mappings is the 
set of bilinear mappings. A function of two variables is said to be bi - 
linear if it is linear with respect to each of the two variables when the 
other is fixed. An example of a bilinear function is given by 

f(x,y) = 3xy . 

This function is linear in each of x and y when the other is fixed 
as can be readily shown: 

f(x,ayj^ + By^) = 3x(ay^ + By^) 

= a3xyj, + 63x72 
= af(x,y^) + 6 f(x,y 2 ) ; 

and 

f(ccxj^ + 6X2, y) = 3 (ax^ + 6x2>y 
= a3x^y + 33x2y 
= af(x^,y) + 6 f(x 2 ,y) . 
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We may note that this satisfies the notion that a bilinear mapping is a 
multilinear mapping with m=2 if both 

+ 6yj^,y2) = at{»(x^,y2) + B’^'(yJ^,y2) (2-2) 

and 

i|)(xj^,ax2 + 6y2) " aif»(x^,X2) + 6^^'(xJ^,y2) (2.3) 

hold for all x^, y^, ct, 3 as noted in the definition of a general multi- 
linear mapping. A bilinear map is sometimes called 2-linear. 

For a given set of vector spaces multi- 

linear mappings from V, x v. x . . . x V to U constitute the set 

12 m 

M(V, :U). The set M(V, ,...,V :U) is a vector space with addition 
1 m 1 m 

defined by 

(^ + 0) (v. , . . . ,v ) ■= 4)(v, ,...,v ) + 6(v^,...,v ) , 

NT ''I’ ’ rn 1 m 1 m 

where the addition on the left side is in M(V, ,...,V :U) and addition 

1 m 

on the right side is in U. Scalar multiplication is defined by 

(v^ , . . . ,v ) = a((j)(v , . . . ,v )) , 

1 m 1 m 

where again (ai) represents scalar multiplication in M(V.,...,V :U) 

1 m 

and ct((j) (v^ , . . . , v^) ) is scalar multiplication in U. Both of these 

principles are sho^>m formally in [11]. Also, it is shown that 

m 

dim M(V, V :U) " n It , 

1 m i 

where 

dim U = n and dim V_, = n^ . 

i i 

Having defined the notion of exactly what comprises a multilinear 
mapping, we are now in a position to consider the idea of a tensor pro- 
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duct. The basic purpose of a tensor product is to convert multilinear 
mappings (in particular, bilinear mappings) into linear mappings. Basic 
ally, this is done because multilinear functions are very complicated ob 
jects and are intrinsically more difficult to handle than linear map- 
pings. We now proceed in defining the notion of a tensor product. We 
consider an arbitrary multilinear mapping 

>^V^x...xV-J-U 
1 z m 

which belongs to M(V, , . . . ,V :U) . It can be shown that there exists an- 

1 m 

other multilinear mapping v e M(V^, . . . ,V^:P) , essentially unique, such 
that there exists a linear mapping p : P -»■ U, which provides 

ij; = y o V . 

The tensor product is said to be composed of P and this multilinear 
mapping v. We now present the formal definition of the tensor product 
[ 10 ]. 

A pair (P,v) is a tensor product of the vector spaces V. , . . . ,V 

1 m 

if the following two conditions are satisfied: 

(1) V e M(Vj^, . . . ,V^:P) and <Im v> = P; 

(2) if U is any vector space over R, and Jp e . . . ,V^:U) is 

arbitrary, then there exists a p e L(P:U) (that is, it is a lin- 
ear map from vector space P to vector space U) such that 

= p o V. 

The property (1) means that all of the vectors generated by v form a 
subspace of P which is equal to P. Property (2) is called the uni- 
versal factorization property , and is expressed by the commutative dia- 
gram in Figure 2.2. It can be shown that the properties (1) and (2) 
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above are equivalent to the following single condition: 

(3) for each multilinear map tfi e . . . ,V^:U) there is a unique 

linear map p e L(P:U) such that 

S3 p o V . 

By the way of shorthand notation, we will define the tensor product ® 

as the multilinear mapping v such that 

v(v, ,v„, . . . ,v_) = V. Q V- a . . . ® V . (2.4) 

1 Z m 1 z m 

Therefore, it follows that the space P may be expressed by a V 2 ® 

... 0 V^. Rewriting Figure 2.2, we obtain the usual form of the commuta- 
tive diagram shown in Figure 2.3. It is shoxTO in [11] that for arbi- 
trary vector spaces V. , . . . ,V a tensor product (P,v) always exists. 

X m 

Also, it can be shown that the tensor product is unique up to an isomor- 
phism. 

We now wish to consider the properties of the space ® ■ 

a V . Namely, the basis and dimension of this space will be examined, 
m 

As was mentioned previously, the bilinear case is nothing more than a 
special case of the multilinear case; so we will first consider construc- 
ting a set of basis vectors for the space 0 V assume that 
and V 2 are both spaces of finite dimension, and that dim = n and 

dim V« = p. We also assume that {e,,e_,...,e } is a set of basis vec- 

2 ^ 1’ 2 n 

tors for V. and that {f,,f f } is a set of basis vectors for V_.. 

X 1 z p z 

It is shown in [1] that the tensor products 

e^^ a f^ , where i = l,2,...,n and j = l,2,...,p , (2.5) 

form a system of linearly independent vectors in the space a 
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Vi X V2 X ... X Q V2 O ... s 



Figure 2.3 Introduction of the Tensor Product Symbol 
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which has dimension n • p, and therefore constitute a basis for the 
space. In this proof, it is noted that 

dim (V^ ® V^) « dim - dim . (2.6) 


If these results are extended to an arbitrary multilinear mapping, then 

we may construct a set of basis vectors as follows. Given that dim 

e and that {e^^^ }, i^^ *= l,2,...,n^ is a basis for k » 1,2,..., 

m, then the space V, © V. ® ... ® V has a basis which consists of 
’ *^12 m 


e, . ® , ® . . . a e, , and 

li, 2i- ki 

12 m 


m 


dim (V, ® V- ® . . . 0 V ) = n dim V. . 
12 m i 


( 2 . 7 ) 


An example will no doubt clarify the use of multiple indices in this 
case. Consider the case where m=» 3 , n^ = dim = 2 , ^2 ^2 = 2 , 

and n^ = dim = 3 . The dimension of ^2. ^ ^2 ^ ^3 2 • 2 • 3 = 12 , 

so there are accordingly 12 basis vectors, which may be listed as follows: 


11 


®21 

® 

^ 31 ’ 

®11 

® 

^21 

0 

® 32 ’ 

®11 

® 

®21 

® 

^ 33 ’ 

®11 

® 

®22 

® 633^, 

11 

® 

®22 

<s 

® 32 ’ 

®11 

® 

®22 

® 

® 33 ’ 

®12 

® 

®21 

® 

® 31 ’ 

®12 

® 

®21 

0 633, 

12 

® 

^21 


® 33 ’ 

®12 

® 

^22 

® 

^ 3 V 

®12 

® 

^22 

0 

® 32 ’ 

®12 

0 

®22 

® ^ 33 > 


For future results and numerical analyses, we shall place much em- 
phasis on the exact ordering of the basis vectors that can be said to 
describe a space. In particular, we shall assume that the basis vectors 
are to be ordered lexicographically . In order to present formally this 
ordering, we will introduce a few elementary concepts from permutation 
group theory [ 12 ]. First, let i^,...,i^ be a set of positive integers 
that satisfy conditions 

1 ^ i- 5 n; 1 5 i„ S n; .... ; 1 S i 5 n . (2.8) 

12 m 
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We shall denote all sequences of these integers by G^. If an integer 

i belongs to the set G°, this integer has m digits, each of which 

belongs to a (n+l)-ary number system excluding zero. If we define 

1g™ 1 as the number of elements in the set g'", clearly Ig*^! “ n"*. 

' n' n ' n' 

The following example is offered to illustrate this concept. Arbitarily, 

3 3 

we choose m = 3 and n = 2. According to our definition, IG2I = 2*=8, 
and the elements may be listed as 

^2 “ { 111 , 112 , 121 , 122 , 211 , 212 , 221 , 222 } . 

This illustrates the concept that the set G^ actually consists of se- 
quences of integers, m integers in each sequence. The range of each 

digit in the sequence is from 1 to n. It can be said that the set G™ 

n 

is ordered lexicographically in the example above, that is, if the ele- 
ments are considered to be an ra digit integer, the sequence of ele- 
ments should start with the smallest number of the base 1) system 
and strictly monotonically increase to the largest number ^ is the 

convention that we will adopt to order our basis vectors . As a final 
example on the calculation and ordering of basis vectors, we extend 
these ideas to multiple tensor products of X and U, which are spaces 
of states and controls. We desire the basis for the space 

X ® X ® U 

where dim X = 3 with basis {x^,X2,X2) and dim U = 2 with basis 
{u^,U2}. The basis is 


(x^ 



0 


^1 

0 


0 

“2’ 

^1 

0 

^2 

0 

“r 

^1 

0 

^2 

0 

^2’ 

^1 

0 

^3 

0 


^1 

<3 

^3 

0 

^2’ 

^2 

0 

^1 

0 


^2 

0 

^1 

0 

^2* 

^2 

0 

^2 

0 


^2 

0 

^2 

0 

^2’ 

X2 


^3 

0 


^2 

0 

^3 

0 

^2’ 

^3 

0 

^1 

0 


^3 

0 

^1 

0 

^2’ 

X3 

0 

X2 

0 
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X3 ® X2 ® U2, X3 0 X3 ® X3 ® X3 ® U2>. 


This is a total of 18 basis vectors. Note that they are listed lexico-' 
graphically. In general, given the repetition p for X and q for 
U, the dimension of 

. X ® ... ® X ® u ® ... ® u 

> — ■ — f \ , I . / 

p q 


is 



(2.9) 


where dim X = n and dim U = m. In this example, note that 

p = 2 , 

q = 1 , 
n = 3 , 
m - 2 ; 


so, 

dim (X ® X ® U) = IG3I • IG2I = (3^) (2) = 18 . 

Thxs is obviously just an extension of the case examined in (2.7), but 
introduces the concept of two different vector spaces which will appear 
in later problem formulations, since we will always be concerned with 
spaces of states and controls. 


2 . 2 Tensor Product o f Linear Mappings f 10 ] 

In this section, we wish to examine the tensor product of two lin- 
ear mappings. These linear mappings will be defined as follows: 

A : ^ , 


B ! ^ Uj . 
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Then, a bilinear map x -»■ can be defined with the action 

i|/(Vi,V2) •=■ (Av^) ® (BV 2 ) 

for € V^, V2 € ^2' relatively easy to verify that ip is bi- 

linear. If we recall the basic definition of bilinearity given in Section 
2 . 1 , for a , 6 e R; w,x e V,; and y,z e V2, 

if)(Vi, ay + Sz) ” (AVj^) ® (B(ay + 0 z)) 

“ (AVj^) ® (Bay + BBz) 

“ (Av^) ® (aBy + BBz) . 

Next, we use the fact that the tensor product itself is a bilinear map- 
ping, so as to obtain 

ay + Bz) = (Av^) ® (BBz) 

“ a(AVj^) ® (By) + B(Av^) ® (Bz) 

*= c'P(v^,y) + B'l^(v^,z) . 

Also, the other half of the proof can be similarly shown, as follows: 
i|j(ocw + Bx, V2) = (A(aw + B^ ) ® (BV2) 

“ (aAw + BAx) ® (BV2) 

«= a(Aw) ® (BV2) + B(Ax) ® (BV2) 

= aij;(w,V2) + B'|/(x,V2) . 

We may express these relationships in the commutative diagram shown in 
Figure 2 . 4 . It is shown in [ 1 ], using the contraction property of ten- 
sors, that the mapping X is indeed a linear mapping, and is equal to 
the tensor product of the two linear mappings A and B. That is, 

X = A ® B , 

where A and B arc as defined previously. This tensor product A ® B 
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is a unique linear map because of the property ( 3 ) of the tensor product, 
vrhich notes the existence of a unique linear map fo each bilinear map, 
given the commutative diagram structure of Figure 2 . 4 . An example will 
clarify these relationships. We shall arbitrarily define the linear maps 
A and B by actions on their basis vectors, and represent these linear 
maps in the usual matrix form. We assume the following sets of basis 
vectors exist for V^, U^, V^, and U2: 

for , Uj^,e2> ; 

for * 

for V2 , {8^,82^ ; 

for U2 , {h^,h2> . 

We next define 

A 6i = - fi + 3 f 2 , 

A e2 = 2 f^ + f2 , 

so that 

Next, we define 


so that 

We now consider basis vectors for and ® U2. Given 

the above set of basis vectors for V^, U^, and U2, we may use 
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0 



Figure 2.4 The Tensor Product of Tt^to Linear Mappings 
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the previous ideas on bases for the tensor product to construct the fol- 
lowing bases: for ® V^, 

(Oj S gj^. 8 82, <=2 8 8^. 02 ® Sj) i 

and for <s U^, 

{I’l ® h^, o h^, f^ ® h^, £2 ® h2> . 

The reader will note that the ordering convention that was previously 
adopted for the basis elements (that is, lexicographic ordering) is em- 
ployed here. We can next construct X “ A ® B by looking at the action 
on the domain basis elements. 

X(Oi ® g^) « (A ® B)(e^ ® 8^) “ (Ae^) ® (Bg^) 

- (~f^ + 3 f 2 ) ® (h^) 

” (-^1 ® + < 3 f 2 ® h^) , 

® So) « (A ® B)(e^ ® g^) “ (Ae^) ® (B82) 

- (-f^ + 3 f 2 ) ® (h^ - 112) 

" (“fj^ ® hj^) + (f^ ® I 12 ) + (3fo ® h^) - (3f2 ® 1^2^ * 


^(©2 ® 8^) “ (A ® B) (g 2 ® “ (Ae2) ® iBg^) 

« (2f^ + f2) ® (h^) 

(2f^ ® h^) + (£2 ® h^) , 

'\(g 2 ® 82) “ (A ® B)(e2 ® 82) “ CAe2) ® (Bg^) 

“ (2fj^ + t'2) ® (h^ - 112) 

« (2f^ s h^) - (2£^ Q I 12 ) + (£2 ® - (fo ® 1^2^ • 

This Implies that the linear mapping X “ A ® B may be represented by 


9 A 


the matrix 
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[X]- 


-1 

0 

3 

0 


-1 

1 

3 

-3 


2 

0 

1 

0 


2 

-2 

1 

-1 


This matrix can also be obtained by the following convention. 


^(Matrix of (Matrix of B) (Matrix of A ) ^2 (Matrix of B) , 

(Matrix of A) 2 j^ (Matrix of B) (Matrix of A ) 22 (Matrix of B) 

This is usually called the Kronecker product of two matrices [13], With 
the above convention, it is possible to verify the result for A ® B, as 


follows . 


A ® B 



-1 -12 2 
0 10-2 
3 3 11’ 

0 -3 0-1 


The Kronecker product will be of much use when deriving the expres- 
sions for the optimal controller in Section IV. In general, it is not nec- 
essary that square matrices be used in computing this product. Given 
(pxq) and (rxs) matrices, the Kronecker product of these two matrices 
is defined to be a (pr x qs) matrix. This can be easily seen as a 
generalization of the above 2x2 case. We also note additional properties 
of the Kronecker product [14]. 


(1) The Kronecker product is associative, that is, (A ® B) ® C = A ® 
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(B Q C) , where A, B, and C are matrices of not necessarily 
equal size. 

(2) (AB) ® (CD) = (A ® C) (B ® D) . This is the factorization property 
of the Kronecker product. 

(3) A ® (B + C) = (A ® B) + (A ® C) . The Kronecker product distributes 
over addition of equal size matrices. 


All of the above properties can be easily verified by using generalized 
matrices and simplifying the result. Note, however , it is not generally 
true that A ® B = B ® A. 


A special case of the Kronecker product occurs when one of the ma- 
trices in the Kronecker product is equal to the identity matrix, denoted 
by 1, , where k specifies the size of the identity matrix. If we assume 

ix 

that A is n X n and B is m x m, then we may define the Kronecker 
sum, A © B, as [13] : 


A©B = A, ®1 +1 ®B. 

m n 


As an example, we consider 



'~3 4 “ 


‘1 2 -6“ 

A = 


, B = 

Oil 


-1 2_ 


3-4 0 


Then, 

A ® 1 = 

m 


fs 4 I 

1 

0 

0 

tl 2j ® 

0 

1 

0 


0 

0 

1 


3 0 0 A 0 0 

0 3 0 0 4 0 ; 

0 0 3 0 0 4 

-1 0 0 2 0 0 

0 -1 0 0 2 0 

0 0 -1 0 0 2 


( 2 . 10 ) 
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1 ® B 

n 


[; 


:] 


2 

1 

-4 



1 2 -6 0 0 O” 

0 1 1 0 0 0 . 

3 -4 0 0 0 0 

0 0 0 1 2 -6 

0 0 0 0 1 1 

0 0 0 3 -4 0 


Therefore, 


A e B = 


4 

0 

3 

-1 

0 

0 


2 

4 

-4 

0 

-1 

0 


■6 

1 

3 

0 

0 

-1 


4 

0 

0 

3 

0 

3 


0 

4 

0 

2 

3 

-4 


0 

0 

4 

-6 

1 

2 


A special case of the Kronecker summation in (2.10) occurs when A 
that is. 


= B, 


A^A = A®1 +1 ®A. 

n n 

We shall find applications of the Kronecker summation idea when we 
consider the derivation of the necessary controller expressions in the 
next section. In order to illustrate partially how the Kronecker summa- 
tion is relevant to system theory problems, we consider the following 
situation [10]. Let us assume a linear dynamical system 

X = Ax 


where A : X ->• X is a linear map and x e X, where X is a vector space 
of finite dimension. Recalling the universal factorization property of 
the tensor product, which allows us to express a bilinear mapping in terms 
of the tensor product and a unique linear mapping, we have the situation 
as depicted in Figure 2.5. The previous results have shown that we can 
consider the linear mapping to analyze the system, since it is known to 
be unique. Since we have already observed that the tensor product is it- 
self a bilinear function. 


27 


(xax) = x®x + x®x 


OF POOR QUALITY 


= (Ax) ® X + X ® (Ax) 

“(A®l)x®x+(1 ®A)x®x 
=[(A®1)+(1 ®A)]x®x 

= (A © A) X ® X . 

Here, the equality x = 1 x was employed, where 1 is the identity map- 
ping of the same dimension as X. One should note the significance of the 
Kronecker sum term, A © A. We shall extend this concept as was applied 
to a linear system to a nonlinear system in terms of the £itate variable x 
and contro.l function u. Since the basic principles involve series ex- 
pansions, va will expect to require higher order terms in the summation. 
The simple example presented above will then be seen as a special case of 
the complete system description which allows for nonlinearities. 


2 . ; Symmetric '* .’ ensor Product 

In this section, we present a brief look at the symmetric tensor al- 
gebra. In order to understand fully this concept in terms of vector spaces, 
it is necessary to first present a few of the ideas concerning quotient 
spaces [10]. Basically, the quotient idea allows us to separate a set 
into two parts: that which is of interest, and that which is not. In 
order to pursue this concept further, we recall the principal ideas of an 
equivalence relation on some set S. Let E be a binary relation on S, 
that is, E is a subset of S x s. We will denote E by the symbol =, 
which is an equivalence relation on S if the following properties hold. 

(i) reflexive, i.e,, s = s for all s e S, 

(ii) symmetric, i.e,, s^ = S 2 implies S 2 = s^, for all s^^. 
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Sj ^ S, 

(iil) transitive, i.e., if 5 S 2 and S 2 = then this implies 
^ ® 3 ’ ® 1 ’ ^ 2 ’ ®3 ^ 

An equivalence relation E divides S into a set of equivalence classes, 
S/ (=) with each class containing elements that are equivalent to each 
other. Information to be discarded is that which would othencise distin- 
guish elements in an equivalence class. If we assume that E is an 
equivalence relation on S, then S/E (read "S modulo E") i-- the set of 
equivalence classes. We define the projection operator, ir as 

7T : S ->■ S/E 


with action 


tt(s) = s’ 

where s e S and s’ is an equivalence class in S/E to which s is 
assigned. If we let f : S T be a function with the property 

s^ E S 2 implies f(s^) = f(s 2 ) > 

then there is a unique function g : S/E T, where T is a set, such that 

g o TT = f . 

These relationships may be expressed via the following commutative dia- 
gram, which is shown in Figure 2.6. This is the "key triangle" that is 
presented in [15]. The main idea in this presentation is that g is 
unique for each f. The existence and uniqueness of this function g 
is shown in [15]. In the algebraic literature, S/E is referred to as 
a quotient set. 


The next step in this sequence is the extension of these concepts 
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on sets to Include vector spaces. Suppose V is an F~vector space, 
and W is a subspace of V, It is possible to define an equivalence 
relation as follows: for "^2 ^ ^ w e W, 

v^ = v^ if v^ => V 2 + w . 

The quotient set of the above discussion is replaced by the quotient 
space V/W, which is also an F-vector space. If we define X as an 
F-vector space and define the linear map 

Q : V -J- X , 

the projection if now becomes a mapping P : V -> V/W. These relationships 
may be expressed in the commutative diagram shown in Figure 2.7. The 
unique linear mapping 

Q : V/W ^ X 

exists if and only if W c ker Q. From Figure 2.7, it is also evident 
that 

Q = Q o P . 


We now can extend these concepts surrounding quotient spaces to the 
subject of symmetric tensors. Namely, we consider the following situa- 
tion. We recall that the tensor product space U ® U has the following 
set of basis vectors, given that a set of basis vectors for U is 




b, ® b. 


9 


0 b^ 


9 


b 


1 


® b 

m 
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b- ® b , 
2 m 


^9 ® » 

2 m 


b ® b 
m m 


(this is a total of m 
basis vectors) 


We next define a linear mapping 


U U ® ... ® U . 
V V / 


(2.11a) 


P P 

This mapping is commonly referred to as the c^nranetrizer [16] and is de- 
fined by 


1 r 

rr I ^ 


(2.11b) 


where a denotes a permutation of variables and the sum is made over all 
possible permutations, with the result being divided by the number of 
permutations. Basically, permutations of indices arise from the inter- 
change of position of these indices. 

First of all, we let Q be a finite set of arbitrary elements. We 

define a permutation on as a one-to-one mapping of onto f?, [12] 

where consists of ra positive integers {l,2,...,m}. We let = 

m 

{k^jk^ , . . . be a set of elements k^ and we let the permutation op- 
erator a be defined as 

^1 2 m 

ct( 1) a (2) a(m) 


o = 
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Here, the first row specifies the domain of the permutation operator a 
and the second row represents ^ image of a. As an example, 


0 


12 3 


|_2 1 3J 

is one such permutation of = {1,2,3} onto itself. Quite obviously, 
there are m! permutations for For the above example, we may list 

all of the members of the image set as 


{123, 132, 213, 231, 312, 321} . 

As expected, there are 3! = 6 members in the set. Returning to the sym- 
metrizer operator, we examine the case where dim U = 2 and p = 3. The 
symmetrizer mapping is 


V ;U®U®U-^-U®U®U. 
s 

The action on the basis elements of U ® U ® U is 

® ® ^1^ “ i" ® b^)) = b^ ® b^ ® bj, , 

(bj^ ®bj^®b-;,+b^0bj^®b2 + bj^®b2®bj^ + b^®b.,®b^ + 
b2 ® bj^ ® b^ H- b2 ® bj^ ® b^) = 

•j (bj^ ® b^ ® b2 + b^ ® b^ ® bj^ + b2 ® b^ ® b^) , 

■^gCbl ® b2 ® b2) = "^ 3(^2 ® ® ^ 2 ^ “ ^^s^^2 ® ^2 ® ^1^ 

■| (b^ ® b2 ® b2 + b^ ® b2 ® b2 + b2 ® b^ ® b2 + b2 ® b^ © b2 + 

^2 ® ^2 ® ^1 ^2 ® ^2 ® ^1^ ~ 

■j (b^ ® b2 0 b2 + b2 ® b^ ® b2 + b2 ® b^ ® b^) , 

^s ^^2 ® ^2 ® ^2^ " i ^^^^2 ® ^2 ® ^2^^ = b2 ® b2 ® b2 . 

For the purposes of this discussion, we will want to consider the 
case where p = 2. We have 
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\ (bj, ® bj^) - b^ 8 . 


POO?? 0 'a>. '7y 


® 1=2^ “ 2 ® bj + b^ ® bj^) . 

n^(b2 ® bj^) = -i (bj^ 8 bj + bj ® bj^) , 

”a<‘’2 ® '=2^ ° ’’2 ® ‘’2 . 

Next, we define the projection ir as follows: 

IT : U ® U •>■ U ® N/ker tt 


The projection operator has the following action on the basis vectors in 


U ® U: 


where 


■rr(bj^ ® b^) = b^ v b^ , 

”<bi ® b^) = b^ V bj , 

IT (b2 ® b^) = b^ V b2 , 

ir(b2 ® b2) = ^2 v b2 , 

bi = ir(b^) . 


The wedge operator v used here is the S 3 nninetric tensor product, which shall 

be defined shortly. If we express the relationships between it and tt^ 

in the form of a coiranutatlve diagram, we have the situation depicted in 

Figure 2.8. Because of our earlier results we may conclude that there 

exists a unique linear map 3 ; U ® U/ker r ^ U ® U with the property 

s 

tr = 3 ° TT . 
s 

If we still assume that the basis elements of the spaces being considered 
here are listed in lexicographic order, then we may obtain the matrix 
representations for ir and it as 
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0 0 

1 1 
2 2 
1 1 
2 2 

0 0 


0 

0 

0 

1 
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and 


[it] = 


1 0 '0 0 
0 110 
0 0 0 1 


Here wa note that the use of the brackets around a linear operator implies 

its matrix representation. Since we choose the basis elements for the 

factor space U ® U/ker it as images of basis elements under tr , then 

s s 

it is possible to determine the matrix representation for 6 as 

10 0 
1 


[g] = 


0 I 0 

0 i 0 


L' 


0 0 1 


It can be easily verified that 


S 


as required. 


A few comments are in order concerning symmetric powers of a vector 


space [16]. If E is an arbitrary vector space, then a vector space 

yP E = E V E V ... V E 




(2.12a) 


together with a symmetric p-linear mapping 
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V : E X . . . X E vP E (2.12b) 

N / 

P 

is called a syrametric power of E if the following conditions are 
satisfied: 

(1) the vectors v ) generate v^E; 

(2) if ip is any symmetric p-linear mapping of E x . . . x E into an 
arbitrary vector space F, then there exists a linear map f : 
v^E -> F such that = f ° v. 

The property (1) means that all of the vectors generated by v form a 
subspace of v^E which is equal to v^E. Property (2) is the universal 
factorization property, and is expressed by the commv.tative diagram shown 
in Figure 2.9. One notes the surprising similarities between the defin- 
ition of the tenso;-: algebra that was presented in Section 2.1 and the 
symmetric censor algebra presented here. As was done previously, condi- 
tions (1) and (2) can be shown to be equivalent to the following single 
condition: 

(3) if is any symmetric mapping of E x ... x e into F, then 
there exists a unique linear mapping f : v^E -> F such that i// = 

f o V. 

It can be shown that the factor space U ® U/ker in Figtire 2.8 is 
isiomorphic to the second symmetric power of U, that is to the space 
U V u. Redrawing Figure 2.8 with this change, we have the commutative 
diagram shown in Figure 2.10. 

We can next define the symmetric tensor product of two vectors s 
and t, with s e U and t e U as 
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s V t = ir(s ® t) . 

As an example, we will assume that dim U = 2, since we have previously 
calculated [it] for this case: 


o 

o 

1 


®1*^1 


'"1*^1 

0 110 


S, t- 


s, t„ + S„t. 



1 2 


12 2 1 

0 0 0 1 


Sjti 


-2"2 


_=2"2_ 



The reader will note that we have expressed 

or dyadic product rearranged as a 4-vector. 

will be explored in the following section. 

find that we may replace a tensor product 

P 4 

u = u © u 0 ... 0 

\ 

p 


the tensor product via the outer 
The reasoning behind this step 
Returning to Figure 2.10, we 

u (2.13) 


with the linear operator 6 acting on the symmetric tensor product, that 


is. 


where 


u 




9 


U = 7T(u) 


and 


~p A 

u*^=.uVuv ... Vu . 

\ / 


As was previously mentioned, this reduction to the symmetric tensor alge- 
bra will be employed whenever we consider the concept of taking deriva- 
tives of the tensor product, a topic which is considered in the following 


section. 
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2.4 Derivatives 

In this section, we wish to examine a few of the concepts that will be 
needed when we derive the results for the various controller terms in Sect- 
ion III. In particular, we will consider the concept of taking a partial 
derivative of a tensor product, which will be needed to solve the funda- 
mental equation of optimality that we are using - the Hamilton- Jacobi- 
Bellman equation. 


Our approach in this section is as follows. First, we present the gen- 
eral definitions of total and partial derivatives and show how these defi- 
nitions can be applied to a very simple case. These ideas are then extended 
to the case in which we are particularly interested, which involves the use 
of the chain rule for abstract derivatives vjhile considering several copies 
of the vector spaces X and U. By the yiay of introduction we may formal- 
ly define the total derivative as follows [17]. If we let V and W be 
normed linear spaces with U open in V, a mapping f : U ^ W is differ- 
entiable at p e U if there exists T e L(V,W) so that for p + x e U, 

X e V, and for | I * | | an appropriate norm on V and W, 

m ||f(p+>t) - , f(p) - ^ 

x^ ||xM 

If such a T exists, then T is unique and called the total derivative 
of f at p, denoted by 


Df(p)x = Tx . (2.15) 

Also, suppose we let V = V, x V , x = (x. ,...,x ) e V and U. 

1 n’ 1 n i 

open in and consider f^ : W with action x^ f (p^, . . . ,p^_^, 

X. , p ,p ) . If f. is differentiable at p, £ U,, we call its 
1 x+1 n X X 1 
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derivative the ith partial derivative of f at p, and denote it by 
D^f(p), p e X ... X u^. Since we will be concerned primarily with the 
problem of taking derivatives of the tensor product, we first examine a 
simple case and note the general concepts involving derivatives. We 
define the mapping f as follows 

f:XxU-^-X0U 

with action f(x,u) = x ® u, where 

X e X and u e U . 

Provided that the mapping f is differentiable, the total derivative of 
f is a linear mapping belonging to L(X ^ U, X ® U) and is defined as 
[17] 

Df (x,u) (Ax, Au) = lim -^[f(x + tAx, u + tAu) - f(x,u)] , (2,16) 

t-K) 


where Au e U and Ax e X are the incremental variables while u e U 
and X £ X are the expansion points. For this particular f mapping 
as defined above, 

Df (x,u) (Ax, Au) = lim^[(x + tAx) ® (u + tAu) - (x ® u) ] 
t-K)J___ 

= lim ~[x ® u + X ® tAu 4- tAx ® u + tAx ® tAu - x ® u] 
t-K) ^ 

= X ® Au + Ax ® u . (2.17) 

Next, we shall be concerned with the idea of partial derivatives. First 

of all, we can define two partial derivatives for this particular case - 

a partial derivative with respect to x and also with respect to u. 

We will denote these by D and D , respectively. 

D f(x,u)(Ax) = lim -^[f(x + tAx, u) - f(x,u)] 

^ t-K) J _ 

= lim — [ (x 4- tAx) ® u - X ® u] 
t-K) ^ 
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= lim -^[x ® u + tAx 0 u - X 0 u] 
t 

t-K) 

= Ax 0 u (2.18) 

Similarly, the partial derivative with respect to u is computed as 

D f (x,u) (Ax,Au) = lim-^[f(x,u + tAu) - f(x,u)] 

t-K) J _ _ 

== lim •— [x 0 (u + tAu) - X 0 u] 

1 » - 

= lim — [x 0 ij + X 0 tAu - x 0 u] 

- t-K) 

=5 X 0 Au . (2.19) 

If we examine the total and partial derivative expressions, it is possible 
to observe the following concepts. Immediately, it is obvious that the 
total derivative is nothing more than the sum of all of the possible par- 
tial derivatives. This is shown formally in [18] for the general case. 
Secondly, we can make the following observations: 

Df (x,u) (0, Au) = D^f(x,u)(au) (2.20) 

and 


Df (x,u) (Ax,0) = D f(x,u)(Ax) (2.21) 

The results that we have presented so far are adequate to study only 
the simple mapping f(x,u) = x ® u. We must next consider the concepts 
necessary in order to take derivatives when there exist multiple copies 
of the spaces X and U, that is, mappings of the form 

f : X X U 0 ... 0 X^0^U 0 ... 0 

V ^ / \« ^ / 

P 9 


In order to examine the general case as stated above, we will need to 
consider certain concepts involving the chain rule in abstract differen- 
tiation [17]. We assume that is open in V^, is open in 

f ; ->■ U^, and g : W for V^, , W Banach spaces. Let p e. 
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be such that f is differentiable at p and g is differentiable at 
f(p). Then g o f ; W is differentiable at p and 

D(g o f)(p) = Dg(f(p)) o Df(p) . (2.22) 

We will assume that h = g ° f. This chain rule for abstract derivatives 
is used in [19] to derive the result that is needed concerning partial 
derivatives of the tensor product. We will assume that X and U are 
vector spaces such that h:XxU-+-(X®X)®U with action h(x,u) = 

X ® X ® u. The mapping h has a derivative at some point p = (x,u) e 
X X U, where x £ X and u £ U. Also, we define the mapping f : X x u 
-> (X ® X) X U with action f(Xru) = (x ® x,u) and the mapping g : (X ® 

X) X u -> (X ® X) ® U with action 

g(x ® x,u) = X ® X ® u . 

Pictorially, we have the commutative diagram displayed in Figure 2.11. 

We also assume that f is differentiable at p and g is differentiable 
at f(p). Then, 

D h(x,u)(Ax) =Ax®x®u+x®Ax®u , 

D^h(x,u)(Au) = X ® X ® Au . 

We will especially be interested in the partial derivative with re- 
spect to the control input, u. If we next allow the mapping h to 
assume the following generalized form, we can extend these ideas pre- 
sented above to include this case. First, we assume that 

h : X X u a- X ® ... ® X a U ® ... ® U . 

N ^ 

p q 

with action 
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h(x,u) =^x ® ... ® X ®^u ® ... ® u, . 

^ ^ V- ' 

p q 

We can compute the required partial derivative as 

D h(x,u)(Au) “X 0 ... ® X. ® (Au 0u® ... ®u. + u®Au ® u ® ... ® u 
u > ^ ^ V ^ / > ^ ^ 

p q-1 q-2 

+ . . . +yU 0 u 0 ... 0 u. 0 Au) 

> — 

q-1 

= x^ 0 D(u'^)(Au) . (2.23) 

This is the result that we will need in Section III to solve the necessary 
minimization problem. Here, as was previously noted, we have assumed 
that X and u are expansion points and Ax and Au are the necessary 
incremental variables. In order to simplify the notational aspects, we 
will drop the "bar" notation when we derive the controller results and 
simply assume that our point of expansion is (x,u). 


Since we have presented some introductory results concerning the 
ideas surrounding the symmetric tensor product, we now consider the par- 
tial derivative operation operating on the symmetric product. We con- 
sider the following situation: 

D^[L(x^ 0 u^)](Au) , 

where L is a linear map operating on the tensor product x^ ® u^. If 
we apply the chain rule for abstract differentiation plus the fact that 
the derivative of a linear map is the linear map itself, the partial 
derivative operator may interchange with the linear map, L, which 
provides 
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“ L x^ 0 D(u^)(Au) . 


u’ = su’ , 


where it is recalled that Q is a linear mapping and 


u =. u V u V ... V u . 
> — 


(2.24) 


Equation (2.24) then becomes 

L[(xP 0 D(Bu'^)(Au)] 

» L[x^ 0 3D(u^)(Au)] 






= L[x^ 0 B(Au V u V ... V u + u V Au v u v ... v u + .. 

\ ^ y \ ^ y 


q-1 


q-2 




+ ^u V ... V u, V Au) ] . 
V 

q-1 


(2.25) 


Next, we will illustrate the specific uses of the symmetric tensor 
algebra structure. It can be shown that all of the terms inside the par- 
entheses in equation (2.25) above are equal. Consider the case where 
q = 2 as an example. Equation (2.25) becomes 

L[x^ 0 B(Au V u + u V Au) ] . (2.26) 

We have previously shown how to compute the symmetric or wedge product 
of two vectors. Using this technique with dim U = 2, we obtain 

/V/ 

AUiUi 


Au V u 


AUjUi + AUj^Uj 
AU2U2 


and 
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Clearly , 


u V Au 




rv r»# 




UoAu, 
2 2 


^fifGINAl 
0^ POOR 


P/IGE ;s 

quality 


rvrf ^ /s/ 

/iu V u = u V Au . 


Therefore, equation (2.25) now becomes 

qL[x^ o 0(u^^ ^ V Au) j , (2.27) 

as there are q identical terms that are added in (2.25) ■ If we recall 


that 


and 


u = 7t(u) 


Au = tt(Au) , 


then (2.27) becomes 


qL[x^ ® 8(ir(u) v ir(u) v ... v ir(u) v 7 t(Au)) 
> ^ 

q-l 


= qL[x^ ® 0ir(u ® u ® ... ® u ® Au) ] 
V — ^ 

q-l 


(2.28) 


= qL[x^ 0 TT^(^u ® u ® ... 0 u^0 Au) ] (2.29) 

S ■ V 

q-l 

= qL[x^ ® TT (u*^ ^ 0 Au) ] . (2.30) 

s 

In obtaining equation (2.29), we note the composition 

TT = 0 o IT 

s 

was used, as was shown in Figure 2.10. Also, we used a property of the 

projection in order to obtain (2.28), namely 

'fr(ai) V TT(a„) v ... v 7r(a ) = Tr(a, ® a„ ® ... ® a ) 

12 n 1 2 n 
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for a. , . . . ,a e U. 

X n 

Basically, we projected into the symmetric tensor algebra in order 
to symmetrize the derivative result, which allowed us to equate all of 
the terms In the summation in equation (2.25). After this operation has 
been accomplished, we effectively return to the (nonsymmetric) tensor al- 
gebra with the Au term in one fixed position. In the following section 
when we begin the derivation for the optimal controller terms, we will 
show that the terms that we wish to differentiate partially are of the 
form presented in (2.24). We will be Interested in showing that certain 
coefficient terms will go to zero for all values of Au, which is why 
these terms must be "factored" out of the expression that is of interest. 

2.5 Discussion 

In this section, we have presented the mathow/; tical preliminaries 
that are necessary to comprehend the remainder of this work. Each of the 
concepts presented so far will be of considerable importance whenever we 
consider the optimal regulation problem in the following sections. Since 
we will be presented with a minimization problem, we will be particularly 
concerned with the procedures of partial differentiation with respect to 
a vector variable. Not surprisingly, the terms that we will be required 
to differentiate will be expressed via the tensor product, hence, the 
reason for presenting the material contained in Section 2.4, which is 
where many of the ideas in the chapter were brought together. 

Thus far, our treatment of the subject matter has been highly theo- 
retical and consistently algebraic in nature, with no reference at all to 
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systems concepts. These concepts are presented in the following section, 
and will appear very similar to some of the classical system theory with 
which the reader is undoubtedly familiar. The difference, of course, is 
the use of the multilinear algebra in order to describe a system, a topic 
that is presented in the first portion of the following section. 



III. A NONLINEAR CONTROL PROBLEM 


The principal idea of this section is to demonstrate how the concepts 
presented in the previous section, on the subject of tensor algebra, can be 
applied to actual control problems involving nonlinear systems. The under- 
lying theme to most of the results presented here is the series expansion 
of a function about a given point. We begin the section with our basic as- 
sumptions concerning the fundamental system description. Having defined 
these preliminaries, we present series expansion ideas as pertaining to 
the eventual solution concepts - that is, in terms of matrices operating 

on basis vectors. A method of calculating the A system matrices via 

pq 

Taylor's series ideas is also presented. The next section focuses on the 
exact procedures that are necessary to derive the needed results, which are 
the controller gain matrices. By construction, this controller is optimal 
in nature, satisfying the Hamilton - Jacobi - Bellman (HJB) equation of op- 
timality. In order to perform the necessary minimization, we shall use the 
concepts presented in Section Ilconcerning the partial differentiation of 
the tensor product with respect to the control variable u. The recursive 
nature of the problem is explored, giving rise to solutions for the con- 
troller terms as well as to terms in the optimal value function. The last 
portion of the section provides a partial verification of the derived re- 
sults as the Linear-Quadratic or LQ problem is verified for the low order 
terms solution. 
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3.1 System Description [1] 


In order to proceed with the system representation in terms of tensor 

til 

expansions, it is necessary to begin with the basic definition of the m 
tensor power of a variable. We define 

u^=u®u®...®u. (3.1) 


J 


V — 

p copies 

The class of nonlinear systems that we will consider are those which can 
be desc“ibed by an ordinary differential equation of the usual form 

x(t) = f(x,u,t) , t £ f^O’^1^ ' (3*2) 

Here, we assume that x(*) ^ is the vector of states and u(0 £ r"* 
is the vector of controls. The systems to be considered here may be 
represented in the following generalized form: 


x(t) = A (t) x^(t) ® u^(t) 

p,'q 

The A terms are linear maps, defined as 
pq 


p+q S 1 


(3.3) 


A : (R*^)^ ® (r"“)^ r"^ 

pq 


(3. A) 


where 


(r“)P a . 


® R 


m 


p copies 


We next define the performance index, J, as 

t. 


J = I M(x(t^)) + i 


L(x(t), u(t), t)dt , 


(3.5) 


where M(x(t^)) and L(x(t) , u(t) , t) are positive convex functionals, 
and L(x,u,t) is assumed to be continuous with respect to t. Also, 


if t^ ->• it is required that the system is asymptotically stable in 
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a sufficiently small neighborhood of the origin. This is the infinite 
time regulation problem, which can be considered as a special case of 
the finite time regulation problem that we will consider. Continuing, 
we assume that the penalty term 

M(x(t^)) = M(x^) , x(t^) = x^ , 

is given by 

M(x^) = I » k S 2 . (3.6) 

The terms are defined as the linear mappings 

: (R^)^ -> R . (3.7) 

We let L(x(t), u(t), t) be given by 

L(x,u,t) = 'l Q (t) x^ 0 u^ , j + k S 2 (3.8) 

j A 

where 

Q., (t) : (R^)^ 0 (R®)^ -> R . (3.9) 

3 ^ 

It is now necessary to define the set of admissible control functions 

as those control functions which can be represented in the usual form 

of a power series in x, that is 

u(x,t) = 'I K.(t) x^ , j = 1,2,3,,,. (3.10) 

j ^ 

where 

Kj(t) : (R^)^ ^ r“ . (3.11) 

The terms follow the previous convention of being linear maps. 

Given the preliminaries presented in this section so far, we are 
prepared to state the basic cptimal control problem, for which a solution 


I* 
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will be constructed in the next section. We desire to find a suitable con- 
trol function u(t) such that for all initial conditions Xq belonging 
to some open neighborhood of the origin of the functional J (per- 

formance index) is minimized. If such a control function is denoted by 
u*(t) and u(t) is any other control function, we require that 

J(Xq,u*) £ J(x^,u) . 

This is to be satisfied for all x^ in an open neighborhood of the origin. 
It is shown in [1] that there exists an open neighborhood of the origin 
such that for all initial conditions in this open neighborhood there exists 
a u(t) that can be represented in the form of a power series as ex- 
pressed in (3.10). Moreover, this solution is unique if the solution to 
the LQ problem is unique. The starting point in the development of this 
control function will be the Hamilton - Jacobi - Bellman equation, or HJB 
equation. It is shown in [20] that the HJB equation is a necessary con- 
dition for optimality. Before presenting these results, however, it is 
first necessary to note a few additional principles surrounding the op- 
timal control problem. We define the so called optimal value function 
V(x(t),t) as follows: 

V(x,t) = min J(x(t) ,u(s) , t) , t_ < t S t^^ , s e [t,tj^] , (3.12) 

uefi 

such that 

X = f (x,u, t) 

is satisfied, with initial condition vector 

xCCq) = Xg . 

This function is called the optimal value function because it is equal to 
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the minimum value of the performance index J on the interval 
Also, the set ^ contains control functions that are expressed as in 
(3.10), that is, 

u(t) = ^ K. (t) x'^(t) , j = 1,2,3,, . . 

3 ^ 

It can be shown that the optimal value function V(x,t) solves the 

following functional equation [20]: 

min V(x,t) + L(x,u,t)] = 0 . (3.13) 

uef2 

This is the functional form of the HJB equation. It also can be shown [1] 
that the optimal value function has the following properties: 

1) If x(tQ) = Xq = 0, then V(x,t) is Identically zero and u(t) = 0; 

2) there is some open neighborhood of the origin in which V(x,t) can 

be represented as a power series in the state variable, x. The first 

term in this series is the quadratic term, so we have 

V(x,t) = I V (t) x^ , k = 2,3,... ; (3.14) 

k ^ 

3) the boundary condition V(x^,tj^) = M(x^) must be satisfied. 

Property 1) listed above is really not surprising if we consider the ini- 
tial conditions as some perturbation from the origin and the state x(t) 

as being an error term. If this perturbation were not present, there would 
be no need for any correction mechanism and the control function would 
be identically zero. Property 2) is very important because it assumes 
that the optimal value function V(x,t) is available as a power series 
in the state variable, x, which is similar to those assumptions made for 
the control function and also for the performance integrand L(x,u,t). 
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The underlying principle surrounding the series expansion of the 
various functions involves the sununation of linear mapping terms oper- 
ating on their respective argument vectors. For example, we assumed 

X “ y A x^ ® u*^ 

P.q ” 

as the given system description. We need to examine the construction 
of the argument vectors x^ ® u*^ which will in turn provide information 
concerning the sizes of the matrices that represent the linear maps 
For the purposes of this example, we assume that dim U = dim X ~ 2. As 
an example argument vector, we consider the tensor product 

X ® u , 

where 


and 


u = [u^ u^] 


X = [x^ x^] . 


In order to form the set of basis vectors for X ® U, we consider the 

T 

outer or dyadic product xu . 

T 



■^r 

II 

<N 

3 

^l“l 

^1^2 


_^2_ 



^2“2_ 


If this 2^2 array is considered to be a 4-dimensional object (that 
is, a 4-vector) with its elements ordered lexicographically, then it 
may be listed as the column vector 

'Vl V2 Vl • 


Alternatively, we can appeal to the bases discussions of the previous 


section. 
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If similar calculations are made for the remaining sets of basis vec- 
tors ^ we can formulate our system description in terms of matrices and 
vectors. Again, we are assuming, for the purpose of this illustration, 
that dim U = dim X = 2, as well as a lexicographic ordering of the basis 


vectors . 


X, 


L^2J L 


X, 


u. 


+ 


10 


^01 


« • 


• • • 


20 


'll 


f" 

2 


mm 


2 ' 

X, 


• • • • 


u. 

1 




1 


+■ 




^1^2 


• • • • 


^"2 

^2^1 


^2'^! 

2 


2 

N) 


“2 





^l“l 


X1U2 


^^2^1 


X2U2 


+ ... 


(3.15) 


Similarly, we can express the series expansions for V(x,t), L(x,u,t), 


u(t), and M(x^) as 


V(x,t) = [ 


V, 


^ 2 

X. 


^1^2 


X2X1 


X, 


+ [ 


V, 


3 

X, 


Xi X2 

^1^2’'l 


X, 


+ . . . , 


(3.16) 


L(x,u,t) = [ 


« • 


• s • J 

‘^1 “ 

+ [ • • • • i 

1- ^ 
“1 

0 

CM 

O' 

^1^2 

CM 

0 

O' 

“1^2 


^2^1 


U2U1 


2 


2 


^2 


_^2 _ 


+ 
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• • « 


‘11 


^l“l 


^l“2 


^ 2'^1 


^2^2 


+ . . . , 


- •" 




•< 


• 

U- 


• • 


X, 


• • • • 

1 

= 



1 

+ 




• • 


^^2_ 


• • • « 


K, 


K. 


X, 


^1^2 


^2^1 


X, 


(3.17) 


(3.18) 


M(x^) = [ 


M, 


1- 2 

^If 

+ [ ] 

r 3 -r 

^If 

^*lf"2f 

“3 

’^If ^2£ 

^2f^lf 


• 

2 


• 

^2f 


3 



^^2f 


+ 


(3.19) 


We should emphasize once again that the representations of these linear 
maps as just illustrated are inherently dependent upon the particular 
ordering convention that is employed for the basis elements. We will 
consistently employ the lexicographic ordering convention in all discus- 
sions. It also would be useful to relate the matrices that we ha\'e 
defined here to the "usual” matrices found in the classical optimal 
control problem. In the traditional approach, a quadratic L(x,u,t) 
would be described as 

L(x,u,t) = x'^Qx + u^^Ru + x"^Cu . (3.20) 

In order to show the relationship, say, between the Q matrix above and 
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the Q 2 Q vector that appeared in (3.17), we assume that dim U = dim X 
= 2, and that Q is represented by 



Then, 


X Qx 


X 2 ] 


^12 

‘^21 ^22 


X, 


Xr 


'^12 ^2 ^21 ^2 ‘‘‘ ^22 ^2 ^ 


‘^12 ^21 ^ 22 ^ 


2 

X, 


^1^2 


X2X1 


X, 


= ^20 " • 


We can similarly determine the components of the other terms as 

needed from their corresponding matrices in the classical optimal control 
problem. 


Thus far, we have defined the A system matrices, but have not 
specified how these matrices may be obtained from a given nonlinear 
system. Basically, there are two methods by which these matrices may 
be obtained. The first is simply obtaining them by inspection. Unfor- 
tunately, this method can only be applied where the variables x and u 
do not appear as arguments for other functions. The following example 
illustrates this method. We assume that we are given the system 

+ 4u^ + 2XjU^ - 5x^2 + , 
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We can calculate the system matrices by inspection. 


» Wi 
• 

X, 


”3 

0” 

* 

X, 


0 

4" 


u. 


0 

1 

2 

0’ 

1 

3S 



1 

+ 



1 

+ 







." 2 . 


1 

0 


_^2_ 


-1 

0 


_^2_ 


0 

0 

0 

0 


^10 ^01 ^11 


x^Ui 


X1U2 


^2^1 


+ 


r* **i 


r 2 - 


r- -1 


r 2 -1 

-5 0 0 0 


^1 


0000 


"1 




+ 




0 ® 
0 
0 
0 

_J 


^1^2 


-2107 


U1U2 

^20 

X2X1 

^02 

“2^1 


2 


2 


to 
I 


c 

to 


For systems where the variables x and u appear as arguments 
for other functions, we must use the other approach, which involves a 
Taylor's series expansion of a function in two vector variables. The 
theory behind this concept is discussed quite extensively in [1] and 
[16] and involves the traditional Taylor's series approach but developed 
in terras of vector-valued tensors and the contraction operator. Another 
approach that defines the series expansion form in terms of tensor pro- 
ducts may be found in [9]. The details of these methods are not really 
relevant to our results here, but the ideas surrounding the definition 
of the various system matrices are quite important. For the purposes of 
this example, we note that we are still assuming that dim U = dim X = 2. 
The system matrices can be computed via the partial derivative op- 

eration, as follows: 
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"10 


'!£i 

3fi‘ 



■«1 

3f,-j 

3x^ 



• *01- 

'“1 

3u2 





3f2 

3Xj^ 


Cx«,u.) 


3uj^ 

3u2j 




'll 


3 ?f. 




2 

3 f, 




3 x^ 3 u^ 


Sx^SUj 


3 x 23 u^ 


3 x 23 u 2 


3 ^f. 


3 ^f, 


3 ^f, 


3 ^f, 


'20 


02 


3Ul 


*X23Ui 

8 X 2802 

* — ! 

CM 

CO 


a^ti 


3Xi^ 

sxj^axj 

SXjSx^ 

3x-^ 

2 

' ^2 

3'f2 

®'*2 


3Xi^ 

Sx^SXj 

3x23x2^ 

8 X 2 ^ 




3^f^’ 








3^f2 

3u/ 

3u^3u2 

3u23u^ 

3U2^ 


^^o’^o- 


(Xo»Uo) 


(Xo,Uo> 


The other values of the system matrices that are needed (depending 

on the problem degree) can be calculated similarly. The point (Xq,Uq) 
is assumed to be the expansion point. An interesting feature regarding 
the above way of calculating the system matrices is how the sequence of 
indices increases lexicographically from left to right in each row. This 
feature is, of course, basis dependent. 

One of the concepts surrounding the HJB equation is the idea of 
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minimization. This minimization will take place after appropriate sub- 
stitutions for V(x,t) and L(x,u,t) are made. Since we are presented 
with a function of two variables, the control u and state variable x, 
we will perform the minimization operation by taking the necessary partial 
derivatives and setting the result to zero. As we shall see, this will 
involve the concepts of Section 2.4 where we examined the possibility of 
taking partial derivatives of tensor products of tXi?o variables. We also 
note that the particular type of control functions that we desire are the 
so called analytic feedback controllers [1]. This requirement is quite 
important in obtaining the final solution. Toward this end, we begin the 
derivatioxi of the controller expressions in the next portion of this 
section. 

3.2 Derivation of Controller Expressions 

la this section, we show how to derive expressions for the controller 
as a function of the system description and the performance index terras. 
The solution will be constructed in a recursive manner by solving the HJB 
equation for the unknown coefficient matrices 

V^(t) , k = 2,3,... 

and 

Kj(t) , j = 1,2,3,... 

such that the boundary condition 

V(x(t^),t^) = M(x(t^)) 

is met. Beginning with the HJB equation 


; 


OI^IGINAL PAGE 13 
OF POOR QUALITY 

we first examine how to compute the time derivative of the optimal value 
function, V(x,t), where we recall from (3.16) that V(x,t) can be ex- 
pressed in sij.ies expansion form as 

2 * 3 

V(x,t) -= V2X + V^x + . . . . 

Therefore, 

V(x,t) = ^2^^ + V2(x ® X + X ® x) + V^x^ + V^Cx ® X 0 X 

+x®x 0 x+x®x®x) + ... (3.21) 

because of the bilinearity of the tensor product and using the ordinary 

product rule. In general, we note the movement of the x then from left 

to right in th> summation and express the general result as 
, “ k 

-3^ V(x,t) = I V, (t) y X ® X 0 . . .yr ® X ® X ® ... ® X 
at ^^2 -“V V 


J-1 

copies 


k-j 

copies 


+ I 


X 


k=2 


(3.22) 


We now examine the substitution for the x terms, namely 
x= y A x^0u*^ , p + qSl. 

^ pq 

p»q ^ 

Making this substitution, and including the expansion tei'-' for L(x,u,t), 
we arrive at the expression in the HJB equation that is to be minimized 
with respect to the vector variable u, namely 

V(x,t) + L(x,u,t) = 


y V, (t) y X «... ® X ® ( y A x^ ® u*^) ® x ® ... ® x + j 

e=2 ^ j=l' p,q Pq SA- k=2 ^ 


X 


j-1 

copies 


k-j 

copies 
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m,n 


(3.23) 


where m+ti ^ 2 and p+q S 1. Recalling the various properties of the 
Kronecker product outlined in Section II we may proceed with the simpli- 


fication. Equation (3.23) then becomes 


I V (t) I I 
k-2 J=’. p,q 


k=2 



® 

1 

‘pq 


X 

k 



X 

+ 

y Q X 



'■ mn 



m,n 




(3.24) 


where we define the multiple Kronecker product of the identity matrix 

1 


1 k-j 

as 


X 


= 1 ® 1 ® . . . ® 1 


■v~ 


3' 


k-j copies 

It would be very advantageous to be able to alter the ordering of 
the terms in the above expression and group the x^ , x and x^ 

terms to obtain a term x^ [16]. However, in order to accomplish 
this, it will be necessary to Introduce a new operator. We define the 
permutation matrix in a slightly modified manner from the definition in 
[21], but with similar results. Let 


where 


S (x ® u ) = X ® u ® X , 

pq 


r p . 


(3.25) 


If we make this substitution in equation (3.24), we have 


dt 


IV , - 

V(x,t) + L(x,u,t) = I V^(t) I I ® ® 1 


k-j X ^k- j 


k=2 


j=i p»q 

k 


pq 


)s 


Iq 


(x ® u^) + y V Tt) x^ + y Q X ® u , 

, k ^ mn 

k=2 m,n 


(3.26) 
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where we define the index Z as 

Jl = p + k - 1 . 

Also, we note that (r=0) is trivially the identity matrix, of dimen- 

Sion equal to the dimension of 

X ® ... ® X 0 U ® ... ® U 




which is equal to 


for dim X = n and dim U = m. 


Z q 
n • m 


An example will be useful in observing the action of this permutation 
matrix. For the purposes of this example, we shall assume that 

k = 2 , 

P = 1 > 

and 

q = 1 . 

Then, we have, from (3.26) 

^2 ^ ^^11 ® V ^21 ^ ® *\ l ^ ■' 

= V2 r(A^« V . 

In general, the operators are linear mappings, and are represented 

in matrix form by a n^ • m'^ square matrix, where dim X = n and dim U 
= m. These matrices have a single "!’* ir. each column and basically re- 
order the listing of the elements in a particular basis vector. For ex- 
1 

ample, can be represented as 
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10 0 0 
0 0 10 
0 10 0 
0 0 0 1 
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For higher order ' p and q, there is an ambiguity that must be resolved 

, 1 

before we can continue. As an example, we consider the operator ^2^, 
with action 


1/2 

82^1^ (x ®u)=X0U®X. 
In matrix form, we have: 


1 

0 

0 

0 

0 

0 

0 

0*^ 


“xiXiUj_“ 


^l“l^l 

0 

0 

1 

0 

0 

0 

0 

0 


X.x, U- 


X-U,X- 

0 

1 

0 

0 

0 

0 

0 

0 


112 

X,X-U, 


112 
X, U„X, 

0 

0 

0 

1 

0 

0 

0 

0 


12 1 
X,X„U„ 


12 1 
X, u x„ 

0 

0 

0 

0 

1 

0 

0 

0 


12 2 
X*x, U„ 

— 

12 2 
x„urx. 

0 

0 

0 

0 

0 

0 

1 

0 


2 1 1 


2 11 
X-U,X„ 

0 

0 

0 

0 

0 

1 

0 

0 


2 12 
X.X u 


2 12 
x„u X 

0 

0 

0 

0 

0 

0 

0 

1 













U 2 2 2- 


L 2 2 2j 


The ambiguity arises in the above example in the construction of rows 2 
and 5 , as well as rows 4 and 7 . We shall adopt the convention thar the 
X vector componenLii will remain in the same order whenever the tensor 
product is reordered. Thus, in row 2 of the above matrix, we choose the 
product equate to x^u^X2» rather than the product X2X^u^. 

This convention will be consistently employed throughout the remainder 
of the work. 


Now that we have clarified the necessary issues, we are ready 
to proceed with the solution. We recall th;^*' the first step in solving 
the optimal regulation problem will be to c i'c'; ate 
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D [^4 +L(x,u,t)] (Au) “ 0 . (3.27) 

u at 

It is necessary to consider what happens whenever equation (3.26) is par- 
tially differentiated with respect to the control variable u. So that 
we may differentiate term-by-term, we shall first list a few of the ini- 
tial terms in (3.26): 


+ L(x.u,t) - + V 2 ((A^g ® + (1^ ® A(,j^)] 


U- ^ ^ XU XU UX A XX 

X ® u + [(Aj_^ ® + (1^ ® A^j)]x2 ® u + ((Ao 2 ® '^x ® ^02>1 

t ® + [(Ajg ® y + (1^ ® A2 o>)x^ +...} + VjX^ + V3{A^q ® A^g ® A^g)x^ + 

on 10 


^*11 ® '^^x ® ^11 ® y®31 ® ® “ + •••> + ••• 

2 2 

+ ^ + Qon X + Qtt X ® U + ... . 


In listing these terms, we have used the Kronecker summation that was pre- 
viously defined in Section II that is, 

A®B = A®1 +1„®B, 

m i 

where A is £ x £ and B is m x m. This sum is well-defined because 
A ® 1^ is an 5,m x Zm matrix as well as 1^^^ ® B. By using this defini- 
tion, we can write, for C n x n, 

A®1 ®1 +1„®B®1 +1„®1 ®C 

m n 5. n 5. m 

= (A®1 +1„®B)®1 +1„®1 ®C 

m £ n £ m 

= (A e B) ® 1 + (1 ® 1 ) ® C 

n S, n 

= A ® B « C . 


This is usually called the multiple Kronecker summation [1] and can be 
viewed as an extension of the normal Kronecker sum. 


If we recall the process of partially differentiating linear map- 
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pings operating on tensor products that was outlined in Section 2 <4, we 
can differentiate (3.28) term by term; 

' 'V ® y ® X ® . 

Du X ® u] (Au) = X ® Au , 

and 


°u *-^02 ^ 

= Qq2 D(u^)(Au) 

“ 2 Qq2 ® • 


If similar calculations are carried out for the remaining terras in (3.28), 
we may list the result as follows; 

Du V(x,t) + L(x,u,t)] (Au) = {2 ® l^) + Q]_^ 

''2 '(^01 ® y ^ ® X ® + (2 Qo 2 "s ® 

+ 3 Q„ (K^ ® ® 1^) + 2 Q^2 [1^ ® x^ (K^ ® 1^)] + + 

V2 ® y + (1^ ® A^^)) + 2 [(Agj ® 1^) sjj + (1^ s A32)] 


'"x ® \ «1 ® "u” “''Ol ® lx'> <^x ® ^01 ® y =21 + 

(1 ^ ® A„. ) ] } x^ ® Au + . . . = 0 
X 01 


(3.29) 


In order to clarify the origins of the terms listed in equation (3.29), 
refer to Figures 3.1 and 3.2, which show the original terms of equation 
(3.28) and their values after differentiation. 


The HJB equation thus far has been solved for the optimal control, 

which can be expressed as a function of x and t. We may write 

u*(x,t) = I K*(t) x^ , j = 1,2,3,... , (3.30) 

j ^ 

where u*(x,t) is the optimal control. If this substitution is made. 
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the resulting equation is (3.29). We will suppress the asterisk notation 
and simply note the optimal gain matrices as (t) . A few comments are 
in order concerning the multiple tensor product of the control vector 
when the substitution (3.30) is made. We will first consider the product 
u ® u. 


u 0 u = (5] K (t) x'O 0 (J IC^(t) x^) 


= I (K (t) 0 K^(t)) X 
Likewise, the general form is therefore 


j+k 


uP = 


3 1 > J 


(K. (t) 0 K. (t) 0 ... O K. (t)) X 




(3.31) 


1 ,^ 2 ^ 


By using these techniques and substitutions, it is possible to obtain 
equation (3.29). 


The factored form of equation (3.29) presents a possible method of 

solution. We desire that this entire expression be identically zero 

for all values of Au. We can express the tensor products x ® Au and 
2 

X 0 Au in matrix form as 


Au , X2 • 


For the case where dim U = dim X = 2, 


^1 = 


x^ 0 

0 X, 


X 2 0 
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X 2 = 


0 




^ 2^1 


X, 


X, 


^ 1^2 


^ 2^1 


X 


2 . J 


and 


Au = 


Au, 


L^'“ 2 J 


Therefore, for an identically zero result, one of the requirements is 
that 

(2 Qo 2 (K, ® 1„) + t(A^^ ® H- (1^ a Ag^)l) x, - 0.(3 

Since the above quantity that multiplies Xj_ is s row vector with four 


components, the general form is 

[Ai Aj ^3 


x^ O' 
0 X. 


X2 0 

0 X, 


= [0 0 ] , 


or 


"a. 

A-' 


X, 


"o' 

1 

3 


1 

“ 

/2 



,^2_ 


0 


which implies that 


.32) 
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Ai «= A2 = A3 = = 0 , 


since A^ is independent of x. This means that our requirement reduces 
to 

2 %2 \ ® '’ll ''2 f^'^Ol ® V®il 'lx ® '^01» = '> • 

It is possible to extend this reasoning to include terms of higher 

2 

order. We can express the tensor product x ® Au in matrix form as 

Xj 4“ • 

This means that we must require 

{2 Qgj (K^ 8 1^) + 3 Qo3 (Kj ® Kj, 9 1^) + 2 Q,, (1^ 9 x, (K^ ® 1^) 1 + 

'’21 1 ''2 '"'ll ® V^21 ^ 'lx ® '^11>> 

2 ’2 "V ® -1 'lx ® -102’' tlx ® ""l ® V’ + 

''3 t'^Ol ® Ix^’sL + 'lx ® ^01 ® + 'Ix^ ® A33)] X2-O ■ ( 3 . 34 ) 

The above quantity that multiplies X 2 is a row vector with eight com- 
ponents, which can be represented as 


f^l "3 ^4 "5 "6 "7 "8^ 


.. 2 


^1^2 


X2X1 


X, 


X, 


^1^2 


^2^1 


^2 J 


= [0 0 ] 


Equivalently, we may write 
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B, B^ B^ B^ 


r 2“] 

X, 


0 

13 5 7 


1 

S3 


h h \ »8. 


^1^2 

^2^1 

2 

Xo 


0_ 



L 2 ^ 




which implies a sufficient condition of 

= 0 , for i 1,2,. ..,8 


In order to show that the above condition is also necessary, we would 
need to answer questions regarding symmetry and the image of a bilinear 
function. Although the components above are Independent of x, we 

cannot argue that it is necessary that they be zero because x ® x can- 
not be made arbitrary. Therefore, we will proceed on the basis of a suf- 
ficiency condition, as proving that the condition is necessary is beyond the 

scope of this work. Since one solution is for = 0, we have 

2 Qq 2 (Kj ® 1„) + 3 Q„3 k, (K^ 8 ® 1^) + 2 (l^- 8 (K^ 8 1^)] + 

«21 + '^2 "^1 ® ^.>=21 ® 2 '' 2 ' t '^02 ® ("x ® '^ 02 >' 


[ 1 ^ 8 8 1 _^)] + V 3 [(Aq 3 8 + ( 1 ^ 8 A 33 8 ys ^3 + 


(3.35) 


It is possible to make similar arguments for the remainder of the coef- 
ficient terms that multiply x™ ® Au, m = 3,4,..., if the tensor product 
is expressed in matrix form as was done for m=l,2 here. For higher order 
cases, it is only possible at the present time to show tne sufficiency 
condition. 


Thus far, we have only considered the partial derivative equation of 
the optimal control u*(x,t). We also must require that [20] 
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_d 

dt 


V(x,t) + L(x,u*(x,t) ,t) 


0 


(3.36) 


along the optimal trajectory. If we collect terms that multiply like 
tensor powers of x in equation (3.36), we obtain the following ex- 


pression: 


+ V2((A^0 e +[(A„, ® ysj, ^ (1^ ® ® K^) + ® + 


'^11 

(1^ ® K^) + Q 

20^=* 

+ {V3 

+ 

V‘ 

® ^x 

’=21 

4 

'^x ® ^01 

® 1 

X 

)S2i + 


(1 ' 
X 

® Aq^) 1 (1^ 

®K^) 

+ (A^j 

1 * 


1® + 

' V 

[(A 

•01 ® lx>=il ® 

(^x® Vi 


® Kp + [(A^^ 

® 1 

X 

)S ^3 + 

(1 

® 

X 

A,i)]a^« 

V 

4 * 

[ (A^„ 0 1 
' 02 X 

:>®12 

4 (1 

X 

0 Aq2)] 


® ® K^) + 

[ (A2Q 

® V 

+ 

("x 

® A2q) ] } 4 

' %2 

[(K^ ® K2) 4 

(Kj 

® K^)] 

+ 

^11 

( 1 ^ ® K2) + Q 

03 


® 

Ki) 

**■ ^12 ^^x 

® 

® 

K3) + Q23 

(1 ^ 
X 

0 K^) 



+ " + . . . = 0 


(3.37) 


Figures 3.3 and 3.4 provide a summary of the origins of the terms that 
comprise equation (3.37). We observe that there are no coefficient terms 
that multiply the first power of x, hence we begin with the second ten- 
sor power. The reader will note that we have collected all of the coef- 
fient terms that multiply like tensor powers of x. In order for the ex- 
pression (3.37) to be identically zero, it is a sufficient condition that 
all of the coefficient terms that multiply like tensor powers of x must 
vanish. The condition has not been shown to be a necessary one because 
of the fact that the product 


X ® X 


does not span the space 

X ® X . 

Again, the proof of this being a necessary condition will not be cousid- 
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Figure 3.3 
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Figure 3.4 


CR'G’NAL PAGE R 

CF POOR qua;.: nr 

ered here due to the scope of the work. This is basically the same quest- 

m 

ion that was found to exist when we collected terms that multiplied x 
® Au, m S 2. At the mcDisnt, it is believed that the use of the symmetric 
tensor algebra will address adequately this issue. Proceeding , then, 
on the sufficiency criterion, this requirement provides the following 
conditions; 



t 

'^2 

4* V ^ 

2 ^10 

^ ) + [(A 

. ® ' 

01 

1 

X 11 


® . 

^01^ 3 (1^ ® K^) } + 





^02 

(K^ ® K^) 

‘’ll ® 

Ki) . 

+ Q 

20 = “ 



(3.38) 

and 













’ 3 - 


{[(Agi e 

1^)S 

X 

21 

0 A 

01 

® 1 )S 

X 

21-^ 


^‘'oi 

)i(i^ 


+ 


^'^10 


^10^ 

}+V2 [(Ag^® 

v^ii 

+ (1 

0 

X 

V’ 

dx® 

Kj) + 


[(A,i 

® V^21 

+ (1 

X ® ^11^ ^ 

(1^ ® 
X 

K^) + 

^^^02 

0 



:®‘"02 

)] 


® K 

^ ® K^) + [ (A 

20 ® ^x^ 

+ (1 

X 

® V 

]}+ Q, 

02 ’ 

(K^ 0 

Kj) + 

(Kj ® 

K,)] t 

<5ii 


® K 2 ) + 

^03 

(K^ ® 

® K^) 

‘’12 

fl ® 

X 

'‘i 

® K^) 

‘’21 

(1^ ® 
X 



+ Q 30 = 0 . (3.39) 


By considering the previous conditions shown in (3.33) and (3.35) 
together with the above requirements, it is possible to observe the re- 
cursive nature of the problem. The objective is, of course, to solve 
for the optimal controller gains K^(t), j = 1,2,3,.... We can solve 
( 3 . 33 ) for Kj^(t) in terms of V 2 (t), and substitute this expression 
into equation (3.38), yielding a differential equation that can be solved 
for ^ 2 (t) with boundary condition ~ ^2’ will 

be the only unknown. After ^ 2 (t) has been obtained, we can easily ob- 
tain the expression for K^(t). A similar procedure can be employed for 
the paired terms of V 2 (t) and K 2 (t), namely, solve equation (3.35) for 
^ 2 (t) in terms of ^^(t) and substitute this expression into ( 3 . 39 ), 
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yielding a differential equation that can be solved for V_(t) with boun- 

3 

dary conditions “ ^ 3 * Th® controller term K^(t) can then eas- 

ily be obtained. This procedure can continue until a sufficient number of 
controller terms are obtained. For the purposes of the example problem 
that will be considered in the following section, we shall assume that we 
are only interested in calculating K^(t) and K 2 (t) although the proce- 
dure could be extended to yield higher order terms if necessary. The re- 
cursive nature of this problem is evident whenever we consider the various 
term dependencies. It can be concluded from the algorithm that the term 
can be determined by knowing only the terms 
K^(t), K 2 (t),...,Kj^_^(t) and V 2 (t), (t) , . . . ,Vj^(t) 


and the system description, which means that the controller term 
can be obtained by knowing only 


K^(t) 


K^(t), K2(t) 


V2(t), V^Ct) , . . . ,Vj^_^j^(t) . 


This result was presented in [1]. In the next section of this chapter we 
show that the LQ problem which has been extensively studied in optimal 
regulation theory is equivalent to the first set of solutions obtained 
with the methods presented here, that is, the set of equations tha!' pro- 
vide V 2 (t) and K^(t). 


3 . 3 The LQ Problem 

In this section, we wish to show that the LQ problem is obtained 
from our results as a special case, namely the case that results from 
the truncation of the system description to linear terms and the perfor- 
mance index to quadratic terms. This result will verify the equation that 
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was derived in the previous section, which is the equation that expresses 
the relationship between K^(t) and V 2 (t) (3.33). 

The normal formulation of the LQ problem has a linear system descrip- 
tion 

x(t) = A(t)x(t) + B(t) u(t) 
and a quadratic L(x,u,t) in the performance index 

1 1 

J = yM(x(t ))+y L(x,u,t)dt , 

0 

where 

L(X;U,t) = x*^Qx + u'^Ru + x^Cu, M(x^) = x"^(t^)Mx(t^) . 

T 

The reader will note the inclusion of the cross term x Cu. Although not 
usually considered in the classical literature, our methods will treat 
this term as being quadratic as well. It can be shown by using the clas- 
sical approach to optimization problems that the unique optimal control 
u(t) can be expressed as 

u(t) = Kj^(t) x(t) , 

where 

K^(t) = - r"^ (B^ V(t) + c'^) . (3.40) 

The matrix V(t) is known to satisfy the differential equation 

* 1 -IT T 1 -IT 

V + V(A - Y B R C ) + (A - “ C R B)V- 

V B R"^ b'^ V + 0 - ^ C r"^ c'^ = 0, V(tj) = M . (3.41) 

If we were to assume that C = 0, which is typical for most practical 
problems, we obtain the usual matrix form of the Riccati equation, which 
is 

V + VA + A^V - VBR“Vv + Q = 0 . (3.42) 
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Our objective in this section will be to demonstrate that the equa- 
tions produced by the expressions previously derived yield the identical 
equations produced by the LQ problem as considered here, namely equations 
(3.33) and (3.38). We begin with our expression relating and V^, 


which is 

''2'V 


and 


If 


this implies that 
Likewise, if 


this implies that 


Also, if 


01 


^1 = 


)] + Qj^j, + 

= dim 

X = 

“’ll 

^12 

_’^21 

^22J 

*"11 

^12 

®21 

"22 


= 0 . (3. A3) 


R = 

'^1 

1 — 

CM 

M 


- 


L 2 

J J 

^02 " ^ 

^1 ^ 

2 


a. 

my 



^12 

C = 



‘^21 

Ut 

^22 


^11 ^‘^ll ^12 ^^21 ^^22^ 


V = 


^21 ^22 
^22 ^24 J 
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^^21 ^22 ^22 • 

It is easily verified that the expression (3.43) produces the following 

set of equations, assuming that the V matrix is symmetric, which 

T 

lows from the quadratic nature of the solution x Vx. We have 


2ri 

=11 

+ 

2r2 

=21 

+ 

2 

=11 

''21 

+ 

2 

^21 

"22 

+ 

c^ = 0 , 


=12 

+ 

=>^2 

=22 

+ 

2 

=11 

''22 

+ 

2 

"21 

"24 

+ 

=3=0, 

2-^2 

=11 

+ 


=21 

+ 

2 

=12 

"21 

+ 

2 

^22 

"22 

+ 

C2 = 0 , 

2^2 

=12 

+ 


=22 

+ 

2 

=12 

"22 

4" 

2 

^22 

"24 

+ 

c, = 0 . 
4 


If the equations produced by the classical solution (3.40) are compared 
with these, one would find that the equations are Identical. 


The other part of the LQ problem involves the verification of the 
Riccati equation (3.41) beginning with equation (3.38), which is 


%2 ® '’ll ® ■'l’ ^ '’20 = ° 


where 



and 



As was previously discussed, the relationship between the Qq 2 ’ ^h’ 
vectors and the classical weighting matrices R, C, and Q in- 
VO Ives the assumption of an ordering convention on the various basis 
elements. Since we are assuming that the basis elements are to be or- 
dered lexicographically, the Q, . vectors can easily be obtained from 
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these classical weighting matrices. 


It can be observed that the expressions (3.41) and (3.38) are very 
similar. The proof that identical equations are produced is not dif- 
ficult but is very tedious when done with the generalized matrices. 
Therefore, for the time being, we will assume that these two expressions 
provide identical equations. When the example problem is considered in 
the following section, we will verify that the derived expression (3.38) 
produces results Identical to the classical Rlccati equation. 

3.4 Discussion 

In this section, we have derived the necessary equations in order to 
calculate the first and second order controller gains and have presented 
the method by which controller gains of any desired order could be calcu- 
lated. Because of the various term dependencies, we showed the recursive 
nature of the solution, namely the alternating solution of the terms in 
the expansion for the optimal value function and the terras in the con- 
troller gain expansion. The LQ problem was partially demonstrated to be 
a special case of the derived algorithm, namely the solution for V 2 (t) 
and K^(t) for generalized matrices from the dimension 2 case. At this 
point, an example is needed to solidify the concepts presented so far 
and demonstrate the method of solution of the recursive expressions. 
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IV. APPLICATION TO SPECIFIC EXAMPLE 


In this section, we present one of the major contributions of this 
work - application of the concepts presented thus far to an example prob- 
lem with a complete set of calculations. The particular example that 
will be considered was analyzed extensively in [9] with particular em- 
phasis on system modeling and model following with a variety of excita- 
tion functions. Because the models generated were proven to closely ap- 
proximate the true solution, these models will be well-suited for our 
purposes. We will present the ideas and techniques for constructing a 
nonlinear control for this particular example based on the expressions 
derived in the previous sections. Basically, we will calculate both a 
first order and a second order control and present an analysis of the 
resulting equations for each case. 

4.1 Problem Requirements and Formulation 

Before we begin the process of specifying an example and choosing 
appropriate weighting matrices that appear in the cost functional, it 
will be advisable to consider the requirements that we wish to meet in 
formulating a meaningful example problem. We will list some of the 
criteria which we have considered in this choice. 

(1) The eigenvalues of A^p. should be in the left half plane 
(that is, they should nave real parts less than zero); the 
idea here is that gas turbine engine models are typically 
stable. 

(2) The choice (x(t),u(t)) 50 is a solution to the differ- 
ential equation. 

(3) The Q^q and > terms are not zero simultaneously. 
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(4) the composite matrix 



is positive definite. 

In order to clarify the controllability issues, we must first assume 
that the time interval is some finite time t^^. ‘The optimal feedback sys- 
tem will turn out to be time varying as long as the control interval is 
finite. This will turn out to be the case even when the system and cost 
functionals are time- invariant, which shall be assumed for the example 
system. The engineering construction of these time-varying functions can 
easily be done on the digital computer using standard integration techni- 
ques. It can be shown [1,22] that if we let t^ -> ®, then we obtain a 
time- invariant controller for a time-invariant system. Basically, control- 
lability is required here in order to ensure that the cost is finite. 

Using a finite interval, there are several examples noted in [22] where 
an opcimal control is obtained for an uncontrollable system. The require- 
ment that the Q^q and 1^2 terms are not zero simultanec .isly excludes 
the trivial case that would produce the optimal control u(t) = 0, al- 
though these terms could be zero individually. 

Now that the example requirements have been formulated, we are ready 
to choose an example system and demonstrate how the gain matrices may be 
calculated for the first and second order cases. The particular example 
that we will choose was effectively studied and analyzed in [9] with par- 
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ticular emphasis on nonhomogeneous model following. We consider the sys- 
tem 


(x,u) 

^2 



2u 

U 2 cosh (x-|^X 2 ) - e 
e sinh (x^^) - e 


1 

u 


1 


sinh (2x^) - 3 sinh (X 2 
2 

cosh ) + sinh (X 2 ) 


) 


> 


and a performance index 


J = -| M (x(t^)) + I 


L(x,u,t) dt 


where we choose the following values for the weighting terms 

Q 2 Q =[2002] , 

Qq2 =[5005] , 

=[ 6006 ] , 

^2 =[1002] , 
=[ 00000000 ]. 



and M, 
J 


We note that the terms were chosen such that the requirement (4) 

is satisfied, which is necessary for a meaningful example problem. Also, 
we are assuming that the penalty term 

M(x(t^)) = M(x^) 

is required to be a convex nonnegative function which means that M(x^) 
is specified in terms of even powers of x^ [1] • 


Since the operators are not directly available, we must use 

the concepts presented in Section 3.1 where these operators were ob- 
served to be related to various partial derivative matrices, namely 


10 




Sx^ 

3X2 


5*2 

3x- 

9x„ 

1 

2 


X = (0,0) 

u = (0,0) 
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and 




3£, 

3f, 

1 

1 

3ui 

3"2 

3f2 


1 \ 

CO 

3Uj 


u = (0,0) 


0 1 

-1 0 


We can easily observe that the requirements (1) and (2) have been satis- 
fied; the origin in X x U defines a solution, since 


f^ (0,0) = 0 , 


for i = 1,2, and has eigenvalues with negative real parts, which 

assures local stability when u is zero. The remainder of the system 
matrices can be similarly calculated as was done for A^q and as 

follows; 






3Xj^3Uj 

axj^au^ 

3X23Uj^ 

3x23u2 

2 

2 

0 

2 

3f2 







®='2®“2j 


( 0 , 0 ) 

( 0 , 0 ) 


-4 0 0 0 

0 0 0 0 


and 


So " 





3f 

3x,^ 

3x^3x2 

3x23x^ 

3x 

CM 

1 CM 

CO 


•Af 2 

3f 


3x, 


3x^3x2 


0 

0 


3x23x^ 


3x, 


02 











3u29u^ 


, 2 
3u2 







_au/ 



9u23u^ 


3u2^ 



0 

0 

0 

0 




-1 

0 

0 

0_ 

• 
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( 0 , 0 ) 

(0,0) 


( 0 , 0 ) 

( 0 , 0 ) 


These methods provide only the nominal or analytical values of the 
system matrices and may be unsatisfactory due to the truncation of 
som« V. the higher order terms. A more accurate method employs the 
least-squares minimization technique using the Singular Value Decompo- 
sition, which is extensively discussed in [9]. Basically, the results 
obtained for the system matrices via least squares minimization provided 
much more accurate responses than those obtained with the standard linear 
approximation and provided a very close match to the "true" solution. 

In order to test and verify this model, various input signals were em- 
ployed with different frequencies and amplitudes. In all of these cases, 
it was found that the model performed remarkably well, and should be well- 
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suited for the generation of the various controller gains which we will 
accomplish in the next section. For the particular example that we have 
chosen, it is discussed in [9] that a degree-2 model (that is, a model 
that is generated by the methods of identification and contains only the 
system matrices A^q, and Aq 2 ) Is adequate to effec- 

tively describe the system as noted above. The identification and ver- 
ification of a degree-3 model involved more effort and programming time 
and was found to only slightly outperform the degree-2 models, and then 
only for a fraction of the test points. Using the least squares identi- 
fication methods, the following matrices were obtained as the degree-2 
model for the example system: 




-2.001 

-3.009 



A -• 





10 

1.006 

1.011 

> 


II 

r-H 

O 

< 

“ 0.002 

0.997“ 

9 


_-1.000 

0.000^ 


^11 " 

“-4. 150 

-0.074 

-0.048 

-0.176“ 

_-0.007 

0.083 

0.008 

0.102_ 


” 0.239 

-0.0725 

-0.0725 

-0.720* 

A„ = 





20 

_-0.323 

-0.064 

-0.064 

0.359. 


"-0.105 

0.0135 

0.0135 

0.012" 

^02 “ 

_-0. 982 

0.0075 

0.0075 

-0.013. 


We observe that the A^^ and A^^^ operators closely approximate the 
nominal values obtained by partial differentiation. By using these 
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values, we can begin to calculate the required terms in the expansion 
for the nonlinear controller. 


4.2 Calculation of Controller Terms 

In this section, we show how to calculate the terms of the nonlinear 
controller expansion. For the sake of simplicity, we will only demon- 
strate the methods to obtain the gain matrices K^(t) and K 2 (t), along 
with the optimal value expansion terms V 2 (t) and V^Ct). As previously 
mentioned, equations for the calculation of higher order terms could be 
developed using the same methods that provided the equations for K^(t) 
and these terms were desired. The first step in the procedure 


^02 


^1 = 


'll 


for 

h 

in 

terms of 

Vj. we 

> 

o 

0 1 

)s 

X 

]■, + (1 ® 
11 X 

V’ = 

;ions , 

that is , assume 

that 

[5 0 

0 5] 




0 

0 

0 ” 



1/2 

1/2 

0 



1/2 

1/2 

0 

9 


0 

0 

1 




^ 2 ' 




> 

^22 




[6 0 

0 6] 

9 
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and use the previously noted values for the system matrix, we find 

that we can express the components of by the following expression: 




’-.0004 

.2 

o' 




- .6 

h2 


-.0004 

0 

.2 


''22 

+ 

0 

hi 


-.1994 

0 

0 


* 

CM 

> 

1 


0 



0 

-.1994 

0 



-.6 

22_ 







U ^ 


In obtaining this equation, we have assumed that = ^ 23 * which cor- 

responds to assuming that the V matrix from classical optimal control 
theory is symmetric. A quick comparison to the equation (3,40) reveals 
that our results are indeed correct. 


The next step in the algorithm is the substitution of equation (4.3.) 

2 

into the equation that was obtained by collecting terms that multiply x 
when the HJB equation is evaluated along the optimal trajectory. This 
result is equation (3.38), which is 


Vj + Vj {(Aj^g ® + ((Agj^ 9 

Qo 2 (K^ ® K^) + Qu (1^ ® K^) + 020 =°- 
If we make this substitution for and use our assumed values for 

the system and performance index matrices, we can obtain the following 
set of coupled differential equations which can be solved for using 

numerical methods. 


90 



-4.044 




+ 3.212 V 22 + *1988 --0008 V 22 + .2 

-3.6072 V23_ - 1.011 V 22 + 1.606 + .1988 

.0004 -.0004 + .2 = 0 , 

-7.2144 V 22 + 2.022 + .1988 V 22 -.0008 

.2 + .2 = 0 . 


POOR QUALITY 


^22 + -2 


0 


V22- 



(4.2) 


As a check on the validity of these equations, we will compare the equa- 
tions listed in (4.2) to those obtained by the Riccati equation, which 

was described in (3.41) as 

’ 1 -IT 1 -ITT 

V + V (A - ^ B R C ) + (A - Y B R C ) V - 

-IT 1 -1 T 

VBR-^B V + Q- ^CR C =0 . 

After sui.stituting the appropriate values .for the known quantities, we 
obtain the following set of differential equations for V(t). 

-4.044 + 3.212 V22 + -1988 -.0008 ^ .2 + .2 = Q , 

V22 -3.6072 -1.011 V22 + 1.606 + .1988 V2^^ -.0004 V22 “ 

.0004 V2^ V2^ + .2 V2^ ^24 " ° 

V2^ -7.2144 V 22 + 2.022 + .1988 -.0008 + .2 + 

.2 = 0 . 


These differential equations are to be solved with the boundary conditions 


''21 

(t^) 

= 1 . 

'22 

(t^) 

= 0 , 

and 



CM 

> 

(ti) 

= 2 . 

Since the differential equations 

presented in (4.2) agree with the classl 


cal Riccati solutions in (4.3), we may conclude that our results are in- 
deed correct. 
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The solution of the first order control problem involves two basic 
steps. In the first step, L e. set of equations (4.3) are solved for V 2 (t) 
and subsequently K^(t) by using numerical integration from the final 
valr.e of V 2 (t). The second part of the procedure involves calculations 
for the control matrix K^(t), which then follows from the known matrix 
values of on the Interval of solution. For this first step of 

the calculations, two approaches may be taken, eith(=r a solution in terms 
of tensor quantities, or a solution taking advantage of the well known 
procedures for Riccati equations. Since we have previously shown that 
our first order results are identical to those produced by the classical 
Riccati equation, it is used in the first order analysis program, which 
is called FIRORDA [23,24]. This program is lipted in Appendix G of this 
report. Numerical values for and V 2 appear also in this appendix, 

denoted by KA and VA. 

In order to test briefly the behavior of the software for the first 
order feedback, the values of the weighting matrices that comprise the 
performance index were changed to five items their assumed values, and 
these results compared to those produced with the example values assumed 
for these matrices. As the Vdiues for these matrices were Increased, we 
found that the solution V 2 (t) of the Riccati equation became more of a 
time varying gain, with greater initial values than those provided by the 
example problem case. This comparison can be made by referring to Figures 
4.1 and 4.2. In obtaining these arrays, we have somewhat arbitrarily 
assumed an interval of t = 0 to t = 5 seconds, with au integration 
stepsize of 0.1 seconds. These results generally agree with the results 
presented in [22]. Note in these two figures that is denoted by VA. 
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Figure 4.1 V2(t) for Example Performance Index 
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A similar procedure is used for the second order analysis of the ex- 
ample problem as was used for the first order analysis. Namely, since the 
values for K^(t) and V^Ct) are known from the first order analysis 
problem, equation (3.34) can be solved for terras of V 2 (t). 

This equation is 

2 %2 % »2 ® '^2 ® y ^21 ® + 

2 Vj [(a^2 ® <ix ® ® ® y> 

''3 ^ ® ^01 ® lx>“21 + <"x ® ^01>1 - ° • 


We note that we have assumed = Qj|^2 “ ^21 ~ ^ purpose of 

the example problem, that is, higher order than quadratic terms in the 
expansion for L(x,u,t) are zero. Since equation (4.4) cannot be solved 
directly for ^ 2 (t) in terms of V^(t), we must multiply the first 
term in the expansion. 


It should be pointed out that the manner in which K 2 is inter- 
twined into (4.4) is one of the most interesting features of the deriva- 
tion in this report, which avoids as much as feasible the use of sym- 
metric algebra and dual spaces. In future work, we hope to return to 
this point. 

Another interesting feature of (4.4) can be seen in more than one 
place, but perhaps especially in the coefficient of on the third 

line. Notice that this coefficient consists of three terms, each in 
essence a rearrangement of A^^ and two copies of 1^. Symmetric al- 
gebra would simplify such sums. 

Returning to the first term in (4.4), we have 
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Multiplying these expressions, we have 
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where denotes the values of rearranged into an 8-vector and G 

A 

is the 8x8 matrix above. Using this notation, we can solve for K_ in 


terms of in equation (4.3) as 

K2 = - fVjKA^, ® ® ^11^1 2 ® 

'"•x ® ® ® lu>l ''3 "^01 ® ^x^>=L + 

»x ® *01 ® ys^i + (1^^ a Aoi)]} G-1 


(4.5) 


The second order solution can be obtained by substituting from equation 
(4.5) into the equation 

Vj + Vj {[(Aoi a + (1^ a Aoi a H- (l/ e Aoi)](l^' a Ki> 

+ (A^o ® A^o a A^o)) + Vj {[(Aq 2 a ys^2 ^^x ® *0l” <^x ® "P ■> 


'(*11 ® (^x ® *ll>l(lx' ® ’^1> ^ '(*02 


V=12 + (^x ® *02>1 


(1^ a a Kp + [(Ajo 


^x> + (^x ® *20” ‘^02 '(‘'l ® '^2> + 


(K2 0 K^) ] + (1^ ® K2) = 0 , 

which is equation (3.39) with the third and higher order terms in the ex- 
pansion for L(x,u,t) equated to zaro. The method of solution of the 
second order control problem is slightly more complicated than the first 
order solution, but employs the same basic steps. 


4.3 Discussion 

In this section, we have examined in terms of an explicit example 
the computational requirements associated with the computation of higher 
order nonlinear feedback control terms based upon the use of tensor de- 
scriptions. Preliminary experience with such computations has been en- 
couraging. Appendix G contains a listing of SECORDA, the program for 
and V^, whose values appear there also, as K2A and V3A. 
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A major purpose of this example study has been to assess the merits 
of using the symmetric tensor algebra and the concepts of dual spaces in 
the calculation of nonlinear feedback tensor gains. To this end, we 
have employed such ideas as infrequently as possible. The results have 
been revealing. For the most part, one can bypass the ideas in question 
and substitute matrix algebra notions such as tra.^spose. Interestingly 
enough, however, there are perhaps two or three instances, one in (4.4), 
where one is strongly inclined to consider the insertion of more techni- 
cality. 

With this observation in mind, we are making a study of the compu- 
tational requirements involved in [1], where symmetry and duals a; both 
used, in order to determine tradeoffs with the current study. 
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V. CONCLUSIONS 


In the three sections preceding, we have carefully examined the compu- 
tational implications of [1], for the case of an example system which had 
received a rather thorough study in a prior grant year. Here we have to 
point out that [1] makes use of technical tensor machinery involving mixed 
tensors on spaces and their duals, the theory of contractions, and the 
formal use of symmetric tensor algebra. The formal purpose of this study 
has been to assess the computational implications of such nonlinear feed- 
back control theory, and in particular to examine the possibility of sup- 
pressing the explicit use of dual spaces and contractions. One benefit, 
for example, of such a suppression would be an avoidance of the distinc- 
tion between tensors and vector-valued tensors. It has been numerically 
clear for a long time, as explained in earlier reports, that the Implicit 
use of symmetry was going to be a certainty. A question, however, was 
whether the explicit use of symmetry in the algebraic derivations would 
really be necessary. 

Having been as careful as we could to carry out the derivation of 
this report while using such concepts as little as possible, we have made 
some important conclusions. First, we had to dip into formal symmetric 
notation for a brief segment of our derivation, in connection with cer- 
tain differentiations. Whether this was an absolute necessity or not 
tends to be outweighed by the fact that a great deal or effort would 
be required in order to circumvent it. We conclude that the formal use 
of symmetry should be pursued. labile this does result in another level 
of equivalence relations, the terminology and symbolism can be induced 
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without a new order of magnitude of difficulty. It tends to be more a 
question of replacing one tensor language with another tensor language. 
Second, we have encountered steps in the equations where duality may be a 
great help in unraveling a complicated twist of symbols. At this point, 
we assess this ability as important, so that we are initiating a careful 
look at the way in which these computations would change if we made use 
of contractions. Should this ti7ork out with the expected benefits, it may 
serve as a very motivating engineering example of the practical importance 
of duality. While the role of duality in general optimization theory has 
been known for many years, there are many common instances in which it can 
be finessed by standard matrix terminology. Such an example could be, 
therefore, very compelling. 

In addition to the work presented in the body of this report, the ap- 
pendices contain a number of items which we also believe to be milestones. 
Appendices B and E deal with ongoing progress in the computer software 
required for tensor modeling and simulation. This developing package is 
the outgrowth of years of work and for the first time makes it possible to 
do examples on a reasonable time scale. In particular, plans are under- 
way to use this new capability to enhance Lli- efforts on tensor model 
identification and order reduction. We expect to make a more complete 
report on this package within the next six months. Appendix C introduces 
a new mathematical viewpoint on the nonlinear feedback control synthesis 
problem. Basically, the idea is to redefine the vector space structures 
for the inputs and outputs of nonlinear systems. For important cases, 
this can be done in such a way that the nonlinear system becomes a linear 
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system with respect to the new spaces. We have this idea under study,- and 
have already shown that it greatly expands the possibilities for defining 
comparison sensitivity functions on nonlinear cases. We envision the use 
of such sensitivity functions to aid in characterizing the quality of 
order reductions and model identifications. In Appendix D we continue 
our investigation into the role of tensors in controller scheduling. We 
expect that the procedures of scheduling may contain one of the keys in 
resolving the tradeoff questions Involving nonlinear model dynamical 
quality versus useful region. Finally, in Appendix F, we present a recent 
compilation of referc'.ces on design by total synthesis, a sequence of 
studies which trace some of their roots back to the first years of this 
grant. 
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ABSTRACT 

Ttie design of nonlinear coatrol systems for gas turbine engines frequently Involves a combi- 
nation of feedforward scheduling and local, dynamic feedback regulation on the desired final 
responses. Scheduling the feedback dynamics, or adding dynamical tuning to the feedforward 
schedules, creates a class of nonlinear dynamical controllers which Is often classical In na- 
ture, as for example the first few terms In a series expansion. Tensor algiibra provides a 
universal setting within which to parameterlre such representations. Moreover, If such mod- 
els are available for the engine Itself, then there exist feedback control theories baaed up- 
on them. In this paper, a model of tensor type Is computed and tested locally on a digital 
simulation of the QCSE gas turbine engine. 

INTRODUCTION 

The use of local, linear dynamical models In control of gas turbine engXi.os has received a 
great deal of attention In the last ten years. While the lion's share of control action for 
such engines tends to be the result of feedforward schedules, the local feedback applied to 
reach desired response points along these schedules Is of great Importance, In particular, 
careful choice of the local controller dynamics can achieve quick, smooch settling, without 
undesirable overshoots in crucial variables, as for example temperatures In the vicinity of 
turbines. 

Such local dynamics are frequently scheduled also, as a function of a smoothly changing phys- 
ical variable, such as a speed. When this Is accomplished, Che local control dynamics be- 
come nonlinear in nature; and key e.xamples can be viewed In terms of vector fields created 
by polynomic functions of state and control, or, more generally, in terms of power series. 
Tensor algebra provides a universal parameteritation within which to represent such schemes. 
Moreover, there exist feedback theories designed to accommodate plant models based upon such 
representations . 

Accordingly, there is Interest In application studies of censor models. In this paper, we 
provide one such study, on a QCSE engine simulator. 

For background, we coi. aider briefly some censor Ideas and Issues associated with nonlinear 
modeling. A short description of the QCSE engine Itself Is given, and then Che application 
Is discussed In detail. 

TENSOR IDEAS 

We begin our discussion with a brief description of the tools to be employed in the nonlinear 
model formulation. Let V and W be real vector spaces and let (®‘^V,at) be a tensor pro- 
duct for r copies of V, where each integer r Is two or greater. For convenience we de- 
fine S'l-V • V and • R. Then by the unique factorization property of the tensor product 
[IJ, for every r-llnear mapping 

U> : - W (1) 

there exists a unique linear mapping 


*This cork was supported by the .National Aeronautics and Space Administration under Grant 
NSG 30A8. 



P/lGf IS 

QUALin 


X : •*'v -* W (2) 

such thac 1 ; • X • 9*^ for the r-llneer napping 9*^ ! V*^ •'^V. If L(V*^;W) denotes Che 

real vector apace of r- linear nappings fron to W, and L(9^V,H) denotes Che real 

vector spaca of linear mappings from «'^V to V, the implication la chat 

LCa'v.tf) - L(V'(W) (3) 


la a vector space iaomorphiso. 

These notions nay be tied Co Che diacuaalon of aba tract derivatives and Che calculus on nom- 
ed vector spaces. As an introduction, equip V and u vith norms and let Z be open in 
V. Suppose chat the mapping t t Z W is dlf ferenciabln at a point p in Z, in the 
usual sense (see, for example, [2,3]). We denote Che derivacJ e of f t Z W at p by 

(Df)(p) t V * W, (A) 

and note that 

Df : 2 - L(V,W); (5) 


that is, Che derivative mapping (A) is a linear mapping, an element of L(V,W). The notion 
extends for higher derivatives, defined in a recursive fashion as 


(D‘f)(p) 


(0(,3'"^f>)(p) 


( 6 ) 


provided the (r-l}ac derlvaclva is differeiclable, since 

D^f(p) < UV.L(V.W)), 
D^f(p) « L(V,L(V,L(V.W))) , 


(7) 


and so on. It can be shown chat there exist isomorphisms 

L(V^;W) ♦ UV.UV.W)) , 


L(V-’jW) * UV,L(V,L(V,W))), 


( 8 ) 


so chat d‘^((p) can be regarded as an r-llnear mapping - W, up to Isomorphism. Me 
suppress this isomorphism and chink of D^f(p) as Just such a mapping. 

It is now straightforward to establish a connection with the censor ideas expressed above. 

The r-llnear mapping (1), for our purposes given by D*^f(p) : - W, can be compoy lu from 

a linear mapping - W and the universal r- linear censor product mapping - a'^V. 

This connection, facilitated by the isomorphisms (3) and (8) , is explored in the section fol- 
lowing for the case of dynamical system rapresencaclon. 

MODEL STRUCTURE 

Suppose thac Che dynamical system which we wish to model is described by the nonlinear ordi- 
nary differencial equation 


X » f(x,u) (9) 

for f : .T » U - .K, where .X and U are normod real vector space3_of states and controls, 
respectively. Using the notation of the preceding section, let (ic,u) be a fixed point in 
Z open in ,X » U, and suppose chat f : ,X « U - .X is of sufficient smoothness on Z. Then, 
formally, 

f(x X, G + u) - ^ ^ (D^f)(x,5)(x.u)^'‘^ (10) 

k-0 

(fc) 

where (x,u)' “ ((x,u),(x,u) (x,u)) k Claes. We note chat the series in (10) could 

be represented by a finite number of terms together with a remainder term in a standard ap- 
plication of Taylor's formula. Indeed, for practical applications, such as the present pa- 
per, a truncation approximation of (10) is considered, Unfortunately, limitations of space 
forbid discussions concerning such issues as existence of solutions to (9) or questions re- 
lated CO the convergence of (10) . 
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Ir « 

W« now oak* uaa of eh« face chat (D f}(x,u) In (JIO) Is a k-llnear napping, vhlch suggoscs 
s Deans of applying censor produce Ideas, Lee « U) ,a^) be a censor produce for k 

copies of X « U. Then we aay tuke ch« unique faccorlzacion 

D'‘f(x,u) - L^()c,u) • (U) 

- - k k 

where L^(x,u) t a (X « U) ■* X Is a linear napping. Now lac che nocaclon (x,u) denoea 

Che k-fold censor produce of (x,u) wlch Icself. Then upon subsclCuClon of (11) Inco (10) 
we have 

f{;rt-x,S+u) - [ ^ L (;,G)(x,u)'‘. (12) 

k-0 

Ic Is shown In [4] chac Che Individual cema of (12) nay be rewrlccen as, for exonple. 

•ji t.j(x,u)(x,u)^ - Ljq(x,u)x«x + Lj^^U,u)x»u + L|jj(x,u)u«u. (13^ 

Ir chls way che fomal expansion (10) becomos 

** " - - 1 ( 

f(x+x,u+u) “II L, ,(x,u)x » u'^ , (14) 

1-0 J-0 

which forms che scruccure for che nonlinear model. 

As alluded Co earlier. In pracclce che series (14) nay be cruncaced In an approxlmaclon of 
(9). The cask In che model building scheme, Chen, Is Co Idcnclfy che psramecers concalned 
in nacrix reprcseiicaclons of che L^<(x,u}, once ordered bases for che spaces In quesclon 
are chosen. For more discussion of the docalls Involved In such an exercise, che reader may 
wish CO consulc (3,61 . 

OCSE ENCrttE 

The Incene of chls aoccion la co supply a brief Incroductlon co NASA's qcSEE ("Qulxla") 

2ulec. £lean, ^horchaul Exporimencal ^ngino— prior to discussing an application of the mtd- 
ellng methodology described above. The QCSE engine Is an .idviinced curbofan designed specif- 
ically for powercd-llf c , shorc-haul aircraft, and combines several Innovative concepts co 
achieve optimal efficiency with quiet, clean operation [7,3], The eight physical quantities 
chosen as state variables for che system Include two fan speeds, four pressures, and two 
temperatures. A digital controller Is incorporated into the overall design (9), and che 
control Inputs are the main burner fuel flow, the fan pitch angle, and che fan nozzle area. 

For che modeling e.xe clsee of chls study, a detailed d'gital slmulaclon developed for che 
QCSE engine [10] Is employed. The primary Input variable co bo manipulated in the digital 
program Is the percentage power demand, PWW. for testing performance over che entire enve- 
lope of operation. Values of Individual Internal variables are extracted and Inserted at 
various locations within Che program. 

.APPLICATION 

Accencion in che following discussion will center around the formulation of a reduced order 
four-state, throe-control analytical model. The engine states chosen are the combustor dis- 
charge pressure (V40S), che core nozzle total pressure (P3GS) , and che rotor dynamics In che 
form of fan speed (NL) , and compressor speed (NH). All three engine control inputs are em- 
ployed, namely, the main burner fuel flow (WFM) , the exhaust nozzle area (AIS) and the fan 
pitch angle (BETAF). 

Appropriate engine operation for the model Idontlf Icacion involves opening the loop by da- 
acclvitlng che controller and Independently Inserting the Individual control Inputs. In this 
strategy, nonilnearlcies of the plant exist— -which might otherwise be less noticeable hod 
the controller bean oresant In the loop. A further explanation of chls strategy as well as 
an alternate one are presented i,' [11), An Important point to note Is that. In Che open 
loop situation, che choice of input control signals Is critical. This Is due to the fact; 
chac the engine Itself has certain physical limits, which In turn have been Incorporated into 
Che simulator. In reality, exceeding these limits could cause severe damage co the engine, 
an example of which is turbine melt down. 

To aid In the selection of input signals, a family of parametric plots have been constructed 
using QCSEE steady state data from Idle (62.51 PtJRX) to maximum power (lOOS PWRX) . Figure 1 
contains an example of one such steady state plot. From those plots a sac of accopcable In- 
put signals can be seloccod. Acceptable state perturbations can bo selected In a similar 
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(fashion. Anocher iraporcanc (aacure of these staady state plots is that thay suggest regions 
of nonlinearity. From Figure 1 it appears that in the locality of 92Z power daciand the eng- 
ine is nonlinear due to the abrupt changes in exhaust nozzle area and fan pitch angle. With 
this in silnd, we shall establish 92Z as the operating point of the present study. Model fora* 
ulatlons at ocher operating points ai ' currently under investigation. 

The following Is an overview of the identification procedure. The QCSEE simulator is cun, 
closed loop, with a 92Z power de^nd for ten_seconds to settle all transients. This pro- 
duces the equilibrium value (x,u), where x la a four-tuple and u is a three-tuple. 

The initial conditions thus generated form Che point of expansion for the secies truncation 
approximation in the model formulation. Vichln Che digital simulation program the controller 
is disconnected by setting the control derivatives to zero. From the steady st.'^te plots a 
point (X,Q) is chosen on the engine oparstlng line at 93Z power demand. The state vari- 
ables are perturbed x from their equilibrium values where 

X “ X - X. 

Furthermore, Che control variables are manipulated so that a cosinusoidal input of amplitude 
3 is inserted, where a is a three- tuple given by the expression 


a « 0 - u. 


The observed states and inputs are sampled over a six second interval; 100 samples are even- 
ly spaced at .06 seconds, and the difference between these values and the corresponding 
equilibrium values, together with the ordered monomials from the censor product terms (see 
[S]) comprise one of two bloclcs of data necessary for the identification. The second bloc 
of data consists of Che state derivative values which are extracted directly from the simu- 
lator at the given sample race. Through use of these data blocks, the parameters contained 
in matrix repressntatioiis of the Lj^,(x,u) can be identified vie a least squares minimiza- 
tion technique. '' 

Using the above procedure, two models have been identified: a second-degree nonl''’iear model, 

and a first-degree linear model. The linear model has been Identified for uie in comparison 
studies. The second-degree approximation keeps second degree t/nsor products which are as- 
sociated with quadratic terms. Accordingly, such a nonlinear Model is expected to outper- 
form the linear model in a region about the point of expansion. 

A simple error comparison criterion is used in testing the performance of the nonlinear mod- 
s'. versus that of the linear model. Let e^ denote the absolute maximum error in the non- 
linear model solution, as compared to the true simulation solution, over the time range of 
simulation for the ith state variable. Similarly, we define for the linear model error. 

Then is the comparison e^ - c^. Thus, if is negative, Che nonlinear model has 

exhibited a smaller maximum absolute error in the ith state, and in chat sense has outper- 
formed Che linear model. Table 1 contains a list of the state variables, their correspond- 
ing QCSEE variable name, their unit of measure, as well as their corresponding state nota- 
tion x^. Samples of the error comparison for various initial condiclc-n, input amplitudes 
and frequencies are presented in Table 2. All input frequencies are in Hertz. 

The error criterion in Table 2 clearly indicates that the nonlinear model outperforms the 
linear model in a region about Che equilibrium point; however, there exists a better method 
for revealing model performance, namely, trajectory comparison. Consequently, a representa- 
tive number of comparative solution plots have been included in Figures 2-10. Figure 2 of- 
fers a simulation of pressure P8GS for a step response, whereas Figures 3— i Illustrate P8GS 
and ML respectively for the frequency set » • (.25, 0., .5). A simulation of P4GS is 
shown in Figure 5 for an excursion away from the typical engine line of operation, and like- 
wise Figure 6 depicts MH. Figures 7-8 illustrate the speeds ML and NH for a 17, amplitude 
control signal, and finally P4GS and ML are seen in Figures 9-10 with initial conditions of 
approximately 1%. 

COMPUTING ENVIRONMENT 

The software package, developed using the extensive capabilities of the IBM and DEC Command 
Procedure Languages and the strengths of FORTRAN and SPEAKEZY, is divided into two segments 
and tailored (,o utilize effectively existent computer hardware. The interactive nonlinear 
model generation segment is implemented on a Time Sharing Option (TSO) of the IBM 370-168 
computer system, where the memory dependent, and highly computational routines of Che pack- 
age can benefit from uva of the virtual memory and floating point hardware. 

Once a nonlinear model is identified on the IBM 370-168, it is transferred to a DEC PDF 11/ 
60, where, in an inceratlve environment, it can be analyzed and compared to a linear model 
and the true solution. This is accomplished through use of Che nonlinear simulation segment 
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of Cha package. In this manner, che couclnea can use boch ehe graphics capablXlcles of a 
Tekcronlx 402S video cermlnal, and a Versaeec electroscaclc prlncer/ploccer for the display 
of data and comparative, trajectories. 

CONCLUSIONS 

This paper has prestinced an application Illustration of tensor modeling to a digital simula- 
tion of the QCSE engine. For plant modeling prior to feedback control, or for representing 
scheduled controllers over an operating line, the tensor algebra offers a universal parame- 
terization which Is helpful In conceptualization and Identification. The case studied In 
this paper offers support to these conclusions. Further work Is In progress. 
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TABLE 1 Variable Ledger for Figures 2-10 

n 

o : True Engine Response 4 : Linear Model Response * : Nonlinear Molel Response 

x. : P4GS (pal) 
Xj: P8CS (psi) 

x^: NL (rpm) 
x^: NH (rpm) 

u, : WFM (lb /hr) 
u^: A18 (?n2) 

u^: 3ETAF (degrees) 


TABLE 2 Comparison Studies 


Initial State Conditions 

Xi Xj X3 x^ 

Input Amplitudes 

“1 “2 

Input Freq. 
Pi Pj ^3 

Error 

£3 =2 *=3 ® 

0.000 

0.000 

0.000 

0.000 

18.91 0.000 -0.111 

0.0 

0.0 

0.0 

-1.040 -0.018 

- 

10.40 -14.00 

0.000 

0.000 

0.000 

0.000 

18.92 0.000 -0.111 

0.3 

0.0 

0.5 

-0.668 -0,012 

- 

7.04 - 8.49 

0.000 

0.000 

0.000 

0.000 

74.29 -21.20 -0.239 

1.9 

0.9 

1.2 

-0.259 -0.004 

- 

4.62 - 8,11 

0.010 

0.001 

0.010 

0.100 

17.00 -2.000 -0.050 

0.0 

0.0 

0.0 

-0.947 -0.017 


7.48 -11.50 

0.010 

0,001 

0.010 

'.'.100 

18.92 0.000 -0.111 

0.0 

0.0 

0.0 

-1.070 -0.019 

. 

10.60 -14.20 

0.010 

0.001 

0.010 

G.lOO 

19,00 -2.000 -0.159 

1.0 

0.8 

0.5 

-0.198 -0.003 

- 

2.49 - 3.16 

0.010 

0.001 

0.010 

0.500 

10.92 0.000 -0.111 

0.3 

0.0 

0.5 

-0.642 -0,012 

- 

6. as - 8.22 

0.010 

0.001 

-0.010 

0.500 

40.00 -5.000 -0.150 

1.5 

l.fi 

1,,3 

-0.354 -0.006 


1.96 - 6.69 

0.010 

-0.001 

0.010 

-0.750 

-18.80 0.000 0.115 

2.n 

0.(i 

1.5 

-0.307 -0.005 

- 

1.71 - 4.89 

-0.010 

-0.001 

-0.010 

-0.750 

-37.46 0.000 0.223 

1,9 

0.0 

*.3 

-0.202 -0.004 

- 

1.59 - 3.42 

-0.010 

0.001 

-0,010 

0.750 

74.29 -21.20 -0.239 

1.9 

1.1 

1.3 

-0.293 -0.003 

- 

5.89 - 6.55 

0.075 

0.001 

20.00 

10.00 

-18.80 0.000 0.000 

1.9 

0,0 

0.0 

-0.343 -0,007 

- 

3.36 - 3.43 

0.110 

O.GOl 

-20.00 

12.00 

55.78 -14.03 -0.198 

1.3 

l.l 

1.4 

-0.447 -0.008 


4.98 - 5.51 

0.120 

0.002 

50.00 

10,00 

37.35 -6.330 -0.159 

2,1 

l.Q 

1.3 

-0.435 -0.008 

- 

5.06 - 4.91 

0.100 

0.002 

75.00 

5.00 

54 00 -12.00 -0.200 

l,'i 

't.O 

1.3 

-0.401 -0.007 


9.80 - 8.80 
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Figure 3 Table 2, Line 7 



Figure 4 Table 2, Line 7 
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Figure 7 Table 2, Line 11 



121 










2^iginal 

■Of POOf? 


IS 


appendix c 

"Nonlinear Multivariable Design by Total Synthesis" 

M.K. Sain 
J.L. Peczkowski 

Proceedings American Control Conference 
Pages 252-260 
June 1982 


123 


original 

OF Poof? 


IS 

Quality 


NONLINEAR MULTIVARIABLE DESIGN BY TOTAL SYNTHESIS 


Michael K. Sain 

Electrical Engineering Department 
University of Notre Dame 
Notre Dame, Indiana 46S56 


Joseph L. Peczkowskl 
Energy Controls Division 
Ihe Bendlx Corporation 
South Bend, Indiana 46620 


Abstract 

In a recei *■ publi.catlon [1], Bristol has pre- 
sented an application theorist's view of process 
control design as it really exists and has chal- 
lenged others to do likewise for areas within their 
own purview. This paper continues Just such an ef- 
fort [2,3,4] by the authors within the domain of 
nonlinear multivariable control of gas turbine eng- 
ines. Under examination is the fundamental notion 
that linear controller descriptions, obtained from 
local actions of nonlinear objects, may be recom- 
bined to produce global nonlinear control action, 
with sufficient integrity to effect closed loop de- 
sign. Total Synthesis refers to a top-down strat- 
egy of Nominal Design and Feedback Synthesis. This 
paper extends the study of the Nominal Design Prob- 
lem (NDP) to nonlinear cases, and presents a new 
case study of robust feedback synthesis for gas 
turbine control design. 

Introduction 

The idea of describing families of curves by 
their tangents has a rich history in mathematics, 
in science, and in engineering. Consider, by way 
of example, the ubiquitous differential equation. 

More generally, the notions of manifold, tangent 
spaces, and geometry are very much a part of modern 
multivariable systems research. 

Not surprisingly, the same notions permeate a 
great deal of control design in various applications. 
Intuitively, one linearizes a nonlinear dynamical 
system at a sequence of points along lines of oper- 
ation considered desirable by the plant manufacturer. 
A suitably rich sequence of points can lead to a 
correspondingly valuable sequence of linear multi- 
variable systems describing local gains and trans- 
ient behavior of the plant along these operating 
lines. From such a sequence of systems one may 
construct a sequence of controllers which effect 
desirable local motions along the lines. Smooth 
global control is then a function of appropriately 
scheduled feedforward and feedback requests, as well 
as scheduling of local controller gains and dynam- 
ics which determine the approach to such requests. 

Bristol [1] believes that experience and in- 
tuition are the crucial requisites for efficacious 
design, and that the best of design-flavored con- 
trol theories can serve only as an introduction to 
the path followed by engineers with experience. Ac- 
cordingly, Bristol also suggests that one should 
seek theories which extend one's intuition and 
which do not presuppose its replacement. 


The local control theory employed by the au- 
thors in this paper was proposed [5] in 1979 by 
Peczkowskl, Sain, and Leake, with Just such a view 
in mind. Conceptually, the method is founded upon 
the idea of a t^ominal Resign problem (NDP), which 
is independent of controller structure and which 
is _ntended as the first step in a top-down design 
procedure. A thorough discussion is given in [6]. 
This paper treats an extension of NDP to the non- 
linear case. Completion of step one in NDP is fol- 
lowed by a second step, called the ^c?edback ^nthe- 
sls problem (FSP) [6]. A case study oi this step 
may be found in [7], which also contains a full 
list of references. An alliance of NDP with FSP 
is called a £otal Synthesis problem (TSP) . The 
case study following in this paper is pari of a 
continuing assault of FSP for the nonlinear case, 
from the view of design practice in gas turMne 
engine control. 

The section following pro' .des mathemati; al 
preliminaries which precede f liscussion of NT-’ for 
the nonlinear plant. Beyond chat, nonlinear nc -i- 
nal design is defined and characterized, after 
which the paper progresses to design of local con- 
trollers for a turbojet case, and the scheduling o; 
these local controls into a global control. 

Mathematical Pre l tnip,i arie3 

In this section, we consider a bijection b : 

S T from a set S onto a set T, with T ad- 
mitting the structure of an F-vector space. As a 
result of the fact that b is bijective, each 
vector t in T can be represented uniquely in 

Che manner b(s) for an s in Sj and each ele- 

ment s in S can be represented uniquely by 

b~^(t) for a C in T. Here, we have denoted the 

inverse of b in Che usual way, b : T -*■ S . 

The commutative group structure (T, 0) on 

our F-vector space T can be used, together with 
b, to induce a commutative group structure (S, Z, 
e) on the set S. The first step in this construc- 
tion is to define the binary operation G : S x S 
S, We do this as follows. Let (s.,s.) € S * S; 
then 

s^ □ S 2 = b b(Sj^) + b(s2>) , 

where the binary operation + in the right member 
is that on T x x to T. Associativity of Che new 
operation can be demonstrated. Indeed, 

(Sj^ G s,) G s^ “ b ^(b » b ^(bs^ + 682) + bs^) 


- b"^((bSj^ + bsj) + bsj) 

- b"^(bSj^ + (bSj + bSj)) 

“ b ^(bSj^ + b 9 b ^(bSj + bs^)) 

» □ (s 2 0 Sj) . 

The unit e can be chosen to be b”^(0), as is 
apparent from the calculation 

s □ b"^( 0 ) - b"^(bs + b 0 b"^( 0 )) 

“ b ^(bs + 0 ) 

“ s . 

For commutativity of the operation, we exhibit the 
steps 

□ S 2 “ b"^(b(s^) + bCSj)) 

- b"^(b(s2) + b(Sj^); 


Finally, for an element s 
additive Inverse d in S 
verify it by 

in S, wj define an 
to be b'^^C-bCs) ) , and 

3 G d “ b ^(bs + b 

9 b"\-b(3))) 

-1 . 


« b ^( 0 ) 



" e , 

as desired. Accordingly, (S, □, b ^(0)) is a 
commutative group. 

Next, we can use the scalar multiplication op- 
eration F < T -*• T on the F-vector space T to 
induce a scalar multiplication F x S S. To do 
this, we define the scalar multiple fs of s by 
f to be 



- b~^(b(fj^s) + b(f 2 s)) 

- (fj^s) □ (f 2 s) . 


IS 

^^ALity 


Next, observe the property 

(fj^f2)s - b‘^((fj^f2)bs) 

- b"^(fj^(f 2 b 8 )) 

- b"^(f^(b 9 b‘^f2bs)) 

- b"^(f^b(f 2 s)) 

- fj^(f 2 s) . 

Finally , 

Is " b ^(lb(s)) 

- b"^(b(s)) 


« s . 

Thus, (S, □, b"^(0)) has been developed in- 
to an F-vector space S. We summarize this fact 
as a theorem. 


Theorem 1. 


Let b : S -*■ T be a bijection onto the F-vec- 
tor space (T, +, 0). Then (S, □, b“^(0)) is al- 
so an F-vector space, with addition 

3 ^ □ S 2 “ b”^(b(s^) + b(s2)) , 

with additive inverse 

b~^(-b(s)) 

for a vector s, and with scalar multiplication 
fs •» b ^(fb(s)) . 


fs = b ^(fb(s)) 


Remark 


for a pair (f,s) in F x x. Notice chat 

f(s^ □ S 2 > “ b ^(fb(Sj^ C S 2 )) 

= b ^(fb 9 b ^(bSj^ + bs 2 )) 

= b ^(fbSj^ + fbs,) 

= b“^(b 9 b'^fbsj^ + b 9 b“"fbs 2 ) 
=■ b"^(b(fsj^) + b(fS 2 )) 

= (fs^) □ (fs 2 ) . 

Moreover, we can also see chat 

(fj_ + f2>s => b"^((fj^ + f2)b(s)) 

= b"^(f^b(s) + f 2 b(s)) 

= b ^(b 9 b ^fj^bs + b 9 b ^f,bs) 


Frequently, S may be given as an F-vector 
space on a commutative group (S, +, 0). Though 
the set S is common to these structures, the bi- 
nary operations + and G are distinct, as wel.l 
as Che units 0 and b“^(0) and the scalar mul- 
tiplications . 

Relative to the Induced space, the bijection 
and its Inverse assume desirable properties. 

Corollar,' 2. 

Regarded as a function 

b : (S, □, b“^(0)) (T, +, 0) , 

the bijection b is a morphism of F-vector spaces, 
as is its inverse 

b"^ : (T, +, 0) (S, G, b"^(0)) . 

Proof : We have only to examine the defining con- 

structions 
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b(Sj^ □ 82) "bob ^(b(Sj^) + b(s2)) 

" b(Sj^) + b( 82 > ; 

b(f3) “ b • b ^(fb(s)) 

■ Eb(s) ; 

moreover, for each in T, we have a unique 

in S such that " b(s^); and so we have 

also 

b + C 2 '' “ b ^(b(Sj^) + b(s2>) 

- b ^ “ b(Sj^ 0 82) 

“ □ 82 

- (b’^(C^)) D (b‘^(t2)) i 

b"^(ft) - b"^(fb(s)) 

« fs 

» fb”^(t) . 


^f^lGlNAL 
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Figure 1. Nonlinear NDP. 

thac NDP is not a model matching problem, in which 
t would also be given and in which only m in 
Figure 1 would appear on a dashed line. 

Now consider a pair (mj^,tj^) of solutions to 

NDP. Characterization of the set ((m ,t )} of 

a a 

all solutions to NDP is severely hindered by the 
fact thac 

p 0 + m,) (p 0 + (p " m 2 ) . 


Remark 

If R is a ring with identity, then all the 
discussions above generalize to R-modules. 

Next, denote by the set of all functions 

from a set R to the F-vector space (S, □, e) . 
Under polncwise conventions, 

(gf 0 g2>(r) - (gj^(r)) D (g2(r)) , 

(fg)(r) - fg(r) , 

R R 

S becomes an F-vector space also [8], as does T 

under Che corresponding operations induced from T. 

The following section defines the NDP for non- 
linear plants, and uses Che properties above Co 
characterize its structure. 

Nonlinear NDP 

The concept of a nonlinear NDP was outlined in 
[3] by the authors for functions on commutative 
groups. Here we extend the idea. Let R, U, and 
Y denote F-vector spaces of requests to the sys- 
tem, controls to the plant, and responses from Che 
plant. It is useful to visualize these, for exam- 
ple, as function spaces, predicated perhaps on time 
sets. Let p : U -*■ Y denote the plant. If the 
feedback action of the controller is well defined, 
then there will be a function m ; R -*■ U gener- 
ating control actions from requests and a function 
c : R -► Y describing plant responses to requests. 
These three functions must then be related by the 
equation 

t = p 0 m , 

which is presented as a commutative diagram in Fig- 
ure 1, The nonlinear Nominal Design problem is to 
find pairs (m,c) in (U^, Y^ such chat the dia- 
gram of Figure 1 commutes. As usual, we point out 


With the ideas of the section preceding, however, 
this difficulty can bo addressed. 

Let (Y, +, 0) denote the F-vector space of 
plant responses, and let (R, +, 0) denote the 
given F-vector space of requests. If p is a bi- 
j action, we can develop on U the F-vector space 
structure (U, □, p“^(0)) of Theorem 1. RelACive 
to this structure, p and p“^ : Y -► U becomCj^ 
isomorphisms of F-vector spaces. Moreover, (U , 

□, e „) and (T^, +, 0) become F-vector spaces, 

with 

e p(r) - p”^(0) 

for all r in R; and 0(r) « 0 in T for all 
r in R. We can Chen define the F-linear map F ; 
uR X ->■ t!^ by setting Che action 

F(ra,t) ” p 0 m - t . 

This leads to the following result. 

Theorem 3 . 

Let p ; U -»■ Y be a bij action onto the F-vec- 
tor space (Y, +, 0). Then a pair (m,t) is a 
solution to Che nonlinear Nominal Design Problem if 
and only if 

(m,t) c Ker F . 

Proof: The assertion is immediate. It may be 

worthwhile, however, to point out Che F-vector pro- 
duct group structure on * T^ defined by the 
operation 

= (m^ I m 2 , + t2> . 
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Remark 

If one wished, he or she could assume vector 
space structure on U, and define a special bi- 
nary operation 0 on T '< T by p(p“^(t.) + p”^ 
(t^)). 

Remark 

An inverse nonlinear NDP, denoted INDF, can 
be defined from the equation 

-1 ► 

m ■ p • t . 


be replaced by + : U * U •* U, 

Suppose next that «'he plant has an internal 
representation 

X f(x,u) , y g(x,u) , 

with appropriate smoothness conditions associated 
with the functions f : X » U X and g : X x U 
Y. Let X be such that 

f(x,u) " 0 , 

and define 


Remark 

Suppose that p were only surjective. It 
follows that p Induces natural equivalence clas- 
ses on Uj and a projection r : U -*• U/= can be 
defined. Then one has the universal factorization 


Ay “ y - g(x,u) , Au " u - u . 

Then p may be assumed to have a local representa- 
tion given by an Impulse response operator 

C + D 


p “ p o ’T 

for a unique p : n/= -*■ Y, which is a bljection. 
The structure (U/=, □, p“^(0)) can be developed, 
and NDP pursued again. Only equivalence classes 
of controls are determined. 

The existence of plant inverses is of the 
first Importance both in theory and in application 
design. In Che next section, we examine briefly 
Che turbojet engine model which will be used for 
our case study. 

Nonlinear Turbojet Model 

A nonlinear model of a simple turbojet engine 
is shown in Figure 2. It is representative, on a 
small scale, of the kind of nonlinear plant with 
which designers of turbine engines and turbine 
controls deal currently in practice. In essence. 

It is a computer simulation, typically constructed 
by engine manufacturers and provided to control 
manufacturers. The nonlinear turbojet model con- 
sists of three integrators, nine nonlinear func- 
tions, including five bi-variant functions, nine 
multipliers and dividers, and nine summing junc- 
tions. The model describes nonlinear dynamical and 
steady state relationships between three inputs: 
fuel flow, Wg, exhaust nozzle area, Aj , and 

turbine vane position, B, and six outputs: eng- 

ine speed, N, turbine temperature, T4, engine 
thrust, FN, tailpipe pressure, P5, and two 
ocher outputs. We propose to Chink of the nonlin- 
ear simulation model as a nonlinear function p 
from a real vector space of control functions of 
time to a real vector space of plant response func- 
tions of time. 

Locally, with appropriate technical assump- 
tions, the nonlinear plant function can be approxi- 
mated by a linear_map P : J -*■ Y, in the neighbor- 
hood of a pair (u,y) in Che relation defined by 
p. I'Then the plant function p is a linear map P, 
the transformation 

P ^(P(u^) + P(u,)) « Uj^ + u^ , 

for the usual vector space structure (U, +, 0). 
Locally, then, the operation Z : U * U - U can 


or by its transform, say P(s), in the usual way. 

The idea is to use these P(s) to determine 
corresponding local descriptions of the parts of 
the controller, and Chen to schedule these parts 
together into a global whole. 

Remark 


In addition to the case study which follows, 
an accompanying paper [9] discusses some additional 
conceptual issues associated with such schedules. 

For Che following case study, three outputs 
have bean selected for control: engine speed, N, 

turbine temperature, T4 , and engine thrust, FN. 
The nonlinear engine model was identified locally 
at five conditions corresponding to 70%, 80%, 90%, 
100%, and 110% speed levels. The engine input is 
given by 


u - (Wj, B) , 

and the selected engine response vector is 
y » (N, T4, FN) . 

By way of illustration, at 100% speed conditions, 
the plant transfer function P(s) and its inverse 
were found to be: 


P(s) 


5.4(.01s+l) 
.13(1,55+1) 
2. 4(. 295+1) 


SS.K.Ols+l) 
-2.7(.50s+l)(,01s+l) 
68.3(.42s+l)(.Qls+l) 


-2704(.01s+U 
:;i6(. 315+1) 
951(. 015+1) 


(,505+1)(. 015+1) 




.18(.23stl)(.0lsfl) 

-.005(-.25*U 

-.OOQl(.74s*l)(.01s*U 


1 . 7 (.aistl)(. 0075 >n 

..08(.01s»l) 

. 00171. 0l5*l)(. 0135+1) 


•.083(.01s+l' 

.015(.01s*l 

-00015(.01s+l, 


!. 0095+1) 
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Figure 2. Nonlinear Turbojet Model. 
Multivariable Design 

A nonlinear multivariable design method, based 
on Total Synthesis ideas, is described and illus- 
trated. Features of the design method include: 

1. an input-output viewpoint; 

2. design for desired response performance; con- 
trol performance, and sensitivity; 

3. a relatively general system structure; 

4. a systematic way to synthesize the nonlinear 
controller. 


Linear Design 


A general linear system structure which com- 
bines TSP ideas with the idea of Comparison Sensi- 
tivity was presented and discussed in [4] and is 
shown in Figure 3. This system structure provides 
coordinated feedforward inputs u^ and loop com- 
mand req”.ists y^^ 


yr ^ 

> . ^ 

LU 




HT 










Figure 3. Robust Controller Structure. 

to a closed loop control system. The feed- 
forward elements coordinate request commands y^; 

the closed loop assures steady state tracking and 
robustness of the outputs, y. The desired over- 
all system response is designated by T. The 
chosen response of the loop is denoted by 

Important controller elements of the struc- 
ture are G, H, and M, which must be designed 
in an acceptable way so as to produce T and Tj^ 

within specifications. It turns out that three 
key equations govern local design: 


M - P"^T 

(I) 

-1 -1 


(II) 

-1 

H - Tj_ (1-Sj^) 

(III) 


Starting point for the design procedure is a 
nonlinear plant model or simulation such as the 
turbojet engine model shown in Figure 2. First, 
it is necessary to establish the desired steady 
state operating conditions of the plant and deter- 
mine available plant inputs and plant outputs. 
Identification of the nonlinear plant along select- 
ed operating lines then can provide local plant 
dynamics and a set of plant transfer function ma- 
trices P(s) relating inputs and outputs. 

Possibilities for control of plant outputs 
using available Inputs can be checked by choosing 
subsets of square matrices of the plant transfer 
function matrix and determining the existence and 
condition of the corresponding plant 'nverse ma- 
trices. Existence of the plant inverse with good 
condition is necessary and vital to obtain reliable, 
independent control of selected outputs with avail- 
able inputs [10-15]. 


Equation (I) is called the synthesis equation. It 
is used to display all admissible responses (T,M) 
and (T ,M^). Equations (II; and (III) are design 


equations for the forward dynamics G and the 
feedback dynamics H, respectively. Note chat all 
control dynamics are defined by selection of M, T, 
T^ and the comparison sensitivity The sen- 
sitivity is defined [16,4] by (1 + PGH)"^. 


These equations provide a basis to design linear 
control systems directly by specifying local re- 
sponse and sensitivity performance. 


Nonlinear Design 

As observed in Che Introduction, the idea of 
describing families of curves by their tangents has 
a rich history in mathematics, in science, and in 
engineering. The method of phase plane portraits 
was already well developed more chan two decades 
ago [17]. In more modern terms, we say today chat 
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scace space descrlpclons, by ordinary differential 
equations, coincide with vector fields on manifolds 
{18]. Solutions on the manifold are curves tan- 
gent to the vectors of the field. 

What about nonlinear design? The fundamental 
notion used in this paper is that linear descrip- 
tions, obtained from local actions of nonlinear ob- 
jects, may be combined to produce nonlinear action, 
with sufficient Integrity to effect closed loop 
control. For example, if a set of local designs 
has given rise to a family {G(s)} of forward dy- 
namics in die loop, then the goal is to link mem- 
bers of Che family together so as to produce a 
nonlinear function g, chat is {G} -*■ g. Now 
each G may be regarded as giving a local approx- 
imation to a part of the vector field. Under re- 
asonable conditions of smoothness, and with enough 
members in the family, a careful linking could in- 
deed lead CO useful g, over regions of interest. 

To extend this notion, one can consider 
choosing, along an operating line, sets of de- 
sired system responses {M,T}, loop responses (M , 
T. }> and sensitivities {S, }, From these perfur-^* 
miSnee choices, sets of controller matrices (M), 
(G), {H}, and {HT} can be generated via equa- 
tions (I), (II), and (III). All of the linear 
sets may be linked and scheduled as a function of 
plant condition to form nonlinear control elements. 
Thus {T} - t, (M) * m, {G} - g, (H) - h, {HT} 
h ® t. 

Desired steady scace operating schedules, 
cre.'isienc control means and protection limits ace 
also needed to provide ocher practical and func- 
tional features for a nonlinear turbojet engine 
control system. These features transform Che 
linear system structure in Figure 3 to the non- 
linear system structure shown below, in Figure 4. 



Figure 4. Nonlinear Control System. 

The structure embodies key rtlacionships of the 
Total Synthesis viewpoint and provides other basic 
features needed for full range, nonlinear control. 
It is used in the design examples which follow. 

Design Examples 

In this section we Illustrate the foregoing 
synthesis ideas by designing a control system for 
Che nonlinear turbojet engine described in Figure 
1 . Recall chat the turbojet has three inputs; 
fuel flow, Wj, nozzle area, Aj , and turbine vane 

angle B; therefore, three outputs; engine speed, 
N, turbine temperature, T4, and engine thrust, FN 
were selected for control. We want to execute de- 
signs CO achieve specific, beneficial output re- 
sponse strategies and show the effect chat sensi- 


tivity specifications have in resisting plant 
parameter variations. Results are Illustrated by 
time response traces for small step commands and 
by full range acceleration and deceleration trans- 
ients . 

Performance Specifications 

Design a multivariable control system for full 
range acceleration and deceleration capability a- 
long the steady state schedules so that complete 
transients are completed in less than three sec- 
onds. Small signal responses of the system are 
desired to produce: 1) fast thrust response; 2) 

smooch, gentle temperature response; and 3) con- 
venient speed response. All should cake place 
without overshoot and without steady state error. 

These requirements translate into the follow- 
ing kinds of refjponse and sensitivity specifica- 
tions : 


1. Tra.'*k output schedules with zero steady scace 
error . 

2. Accelerate or decelerate from 70% to 100% speed 
levels in less chan 3 seconds. 

3. Local System Responses (T) - Decoupled 

Speed - .5 second lag (? 70% speed 
- .2 second lag 0 100% speed 
Temp. - .5 second lag - constant 

Thrust - .2 second lag - constant. 

4. Local Closed Loop Response (Tj^) - Decouple*'. 

.2 second lag constant for all outputs. 

5. Local Sensitivity (Sj^) 

a) Unity feedback: Sj^ » (I-T^) ; 

b) Ten times better than unity feedback 

The following response and sensitivity matrices at 
100% speed condition were obtained: 


T 

T 

S 

S 


L 

m 

L 

a 

L 


diag 

dlag 

diag 

dlag 


, 

\2s+l‘ 

( 1 , 
\2s+l’ 

'5e 

• - r J — . 
\2s+l’ 

f .02s 

\02s+l’ 


1 

1 . 

.5s+l’ 

.2s+l' 

1 

1 . 

.2s+l’ 

.2s+l' 

,2s 

.2s . 


.2s+l’ .2s+l^ 
.02s .02s 


.02s+l’ .02s+l 


(H-I) 

) 


The forward control dynamics G(s) af 100% 
speed condition with unity feedback are: 


G(s) 


.S0(.23s+1) 

-.025(-,2s+l) 

.0005(.74s+l) 


8,3(.01s+l) 

-.40(.01s+l) 

.0085{.013s+l) 


-.42(.01S+l) 

.075(.0ls+l) 

.0008(.01s+l) 


s{.01s+l) 


Contr.. .ipr elements M, G, H and HT were 
calculated at five engine speed conditions. These 
sets were scheduled as a function of engine speed 








Co form nonlinear control system elements m, g, h 
and h • c. For example, the form of the nonlinear 
controller g so constructed Is shovm In Figure 5 
below. 



Figure 5. Nonlinear g. 

Simulation Results 

Small step transients of the nominal engine 
with sensicivlty feedback system are shown In Fig- 
ure 6, The output responses verify desired small 
signal performance: thrust response is fast (.2 

second lag); temperature response Is smoc^Ch (.5 
second lag), and speed response varies from .5 sec- 
ond lag at 70X Co .2 second lag at 100% speed con- 
dition. Corresponding input responses are shown 
in Figure 7. 

Full range acceleration and deceleration tran- 
sients of Che nominal engine with sensitivity feed- 
back system are shown in Figure 8. The outputs 
track Che requests without overshoot and Che trans- 
ient time is less than three seconds. Correspond- 
ing input responses are shown in Figure 9. All 
Inputs are within desired limits. 

To show the effect of sensitivity specifica- 
tion on plant parameter variations, Che time cons- 
tants of the engine affecting speed, temperature 



Figure 6. Outputs. 

Nominal Engine; Sensitivity F/B. 


INPUTS 

NOM. ENGINE SENSITIVITY F/S 



Figure 7, Inputs. 

Nominal Engine; Sensitivity F/B. 

and tall pipe pressure response were all doubled. 
This produced a nonnominal engine. 

Small step responses of a nominal and non- 
nominal engine with unity feedback control are pic- 
tured in Figure 10. Deviations from desired output 
responses are caused by the engine parameter vari- 
ations . 

Imposing a sensitivity specification which is 
effectively ten times better than the unity feed- 
back design produced feedback dynamics H = diag 
((.2s+l)/(.02s+i)) . Step responses for the sensi- 
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OUTPUTS 



OUTPUTS 

UNITY F/B 



Figure 8. Oucputs. 

Nominal Engine j Sensitivity F/8* 

INPUTS 

NOM. ENQINE SENSIrlVITY F/B 



Figure 10. Outputs; Unity F/B, 
Nominal and Nonnomlnal Engines. 


OUTPUTS 
SENSITIVITY F/D 



Figure 9. Inputs. 

Nominal Engine; Sensitivity F/B. 

tivity feedback design with both the nominal and 
nonnomlnal engines are shown in Figure 11. Devia- 
tions due to engine parameter variations are vir- 
tually eliminated. 

Full range acceleration and deceleration tran- 
sients of the sensitivity feedback controller with 
both nominal and nonnomlnal engines are shown in 
Figure 12. The sensitivity design feedback con- 
troller maintains full range output transients es- 
sentially at nominal conditions, successfully hand- 
ling plant parameter variations. 


Figure 11. Outputs; Sensitivity F/B. 

Nominal and Nonnomlnal Engines. 

Summary Remarks 

A nonlinear control synthesis method based on 
TSP viewpoint was discussed and illustrated. A 
three input/three output turbojet engine example 
demonstrated a feasibility to achieve desired sys- 
tem response and sensitivity specifications. 

A concept of the nonlinear Nominal Design Prob- 
lem (NDP) was presented and discussed, extending 
and building on earlier Total Synthesis Problem 
(TSP) theory and ideas. Additive structure was ob- 
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OUTPUT3 



Figure 12. Outputs; Sensitivity F/B. 

Nominal and Nonnominal Engines. 

talned by a process of in ; icing a special binary 
operation on the control input space. Though not 
a now mathematical idea [18], this concept seems to 
fit constructively into current design developments 
in nonlinear control. 

Research to develop nonlinear control synthe- 
sis methods is needed. It is felt that the input- 
output TSP viewpoint offers possibilities to devel- 
op useful, systematic and straightforward methods 
for nonlinear multivariable control synthesis. 
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CONTROLLER SCHEDULING: A tOSSIBLE ALGEBRAIC VIEWPOINT 


Michael Sain and Stephen Yurkovlch 
Department of Electrical Engineering 
University of Notre Dame 
Notre Dame, Indiana ^6556 


Abstract 


In the applications, one common way to design 
a control system for a nonlinear plant is to local- 
ize its behavior along lines of operation specified 
by the plant manufacturer, to develop linear multl- 
variaole controls for these localizations, and to 
schedule those controls with key plant variables 
which vary smoothly along operating lines. An im- 
portant part of practical design lore, the art of 
controller scheduling has received little modern 
attention from the conceptual point of view. This 
paper describes four basic types of scheduling ques- 
tions and outlines some of the theoretical issues 
associated with them. Schedules are considered in 
terms of state equations; however, some relations 
with the Input/output description are discussed, 
together with an analysis of the effects on the 
overall configuration of approximations made to 
the Individual subsystems. 

Introduction 


Bristol [1,2] has likened the process of con- 
trol design to the use of idioms In a language. 

At least three types of id^jms can be identified. 
First, there are idioms which have been with man- 
kind for such a length of time that they seem uni- 
versal to the human psyche. In some sense, feed- 
back itself is an example of such an idiom, inas- 
much as it may be traced at least back to ancient 
Arabian water clocks. Second, there are idioms 
which are the characteristic of certain authors. 
Several classic examples are the Nichols chart, 
the Bode plot, the Evans loci, and the Nyquist 
plot. And third, there are idioms which are typi- 
cal of certain types of control applications. An 
example is chat of gas turbine control systems [3]. 

Because of the idioms of type three, any ap- 
plication of control design has idiomatic features. 
In a sense, the Cask of Che control designer is to 
blend the idioms of the application with universal 
idioms, with idioms of classical and modern authors, 
and with his or her own idioms, so as to produce a 
melodious and effective composition. 

It goes without saying that some idioms do not 
play well together. In some areas of application, 
this may account for the famous theory/appllcacion 

gap- 

One universal -Jiom is to attack the overall 
system design by breaking it down into manageable 
pieces. An important case of this type of thinking 
arises in the design of certain classes of nonlinear 
systems. Examples in point may be found in the 
area of gas turbine control. For discussion of 


s.vme of Che ideas involved, as well as additional 
references, see the companion paper [4]. In brief, 
the nonlinear engine is linearized locally along 
lines of operation agreed upon by the manufacturer 
and the control contractor. These linear multi- 
variable localizations are used to develop a family 
of local controllers, which are then sewn together 
by scheduling control gains and dynamics with some 
engine variable, as for example speed, which varies 
smoothly along operating lines. 

As pointed out by Bristol [1], Che idioms have 
to blend together. In the case of scheduling, the 
methods used for design of the local, linear mul- 
tivariable controllers have to be amenable to a 
common thread of smooch scheduling, else a global 
whole is not obtained, but only a sum of parts. 

The goal of this paper is to examine in an in- 
troductory way certain of the conceptual questions 
associated with scheduling. In particular, we 
would like to determine something about the struc- 
tures of common scheduled systems, their approxi- 
mations, and how they are affected by interconnec- 
tions one with the ocher. 

What follows should be regarded as exploratory 
in nature. Though we will do some things in con- 
siderable detail, it nonetheless remains true chat 
we will be answering only a portion of the questions 
which we raise. 


Simple Example 

Consider Che elementary dynamical system 

X “ -ax + bu . (1) 

The transfer function associated with (1) is of 
course 


b 

s+a 


( 2 ) 


Rewritten in terms of gain and txme constant, (2) 
becomes 


k 

TS+1 ’ 


(3) 


where 


k « b/a , T » 1/a . (4) 

Suppose chat we wanted to schedule the gain k as 
a function of the input u, say 
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2 

k(u) ■ + Bj^u + Yj^u . 

(5) 

Then the scheduled system would look like 



2 3 

X ■ -ax + aOj^u + aBj^u + ay^^u . 

(6) 

Alternatively, we might schedule Che time constant 
T as such a function, for example 


2 

t(u) " «2 ^2“ ^ ^7“ ’ 

(7) 

In 

which case we would have 



2 

a ■ 1/ + SjU + YjU ) 



-1 -2 
- a, - 62“2 “ + • • • 

(8) 

so 

that 



. -1 -2 

X “ -02 X + 6202 ux + bu + ... 

(9) 

on out to a denumerably infinite number of terms. 
Next suppose that we wanted to schedule the gain 
k or the time constant t as a function not of 
u but of X, in the manner 


2 

k(x) - , 

(10) 


2 

t(x) - + Y4X . 

(11) 


Then Che scheduled systems would be 

X “ -ax + aa^u aS^xu + ay^^'u, (12) 

X “ -a.^x + x^ + bu + . . . , (13) 

4 4 4 

again with a denumerably infinite number of terms. 

Generally ipeaking, the polynomic scheduling 
concept cendt to convert the system (1) into a sys- 
tem of the form 


K " I I c x^u^ . (14) 

i-0 j-0 ^ 

Indeed, if the original system (1) were of Che more 
general form 

i- 1 J a x'^u® . (15) 

k-0 m-0 


and if Che parameters w?7e scheduled in an analog- 
ous way, such as 


km 


I I 

p"0 q=0 


\mi 


pq 


P 

X u 


(16) 


Chen (15) becomes 




I I I I 

k»0 m=0 p=0 q*»0 


ipq 


^k+PuHtfq ^ 


which can be formally rearranged in the same form 
as (15). In broad terms, then, (15) is closed 


under formal power series scheduling. 

Because of this closure feature, we find In- 
terest in systems of this type. The next section 
given a brief, multivariable motivation. 

Multivariable Motivation 

When a nonlinear plant has been linearized in 
Che neighborhood of a point on Che desired operat- 
ing line, the resulting, local, linear multivari- 
able control problem must be resolved by some 
chosen procedure. Moreover, associated with this 
chosen method will be various related and compat- 
ible theoretical viewpoints. It is not the pur- 
pose of this paper to argue the merits of one or 
more of these theories. Instead, we wish to select 
a way of thinking which is comprehensive enough to 
embrace the thoughts of niimerous approaches to such 
conceptualization. In this way, the scheduling 
classes which we introduce will, hopefully, be 
broad and representative of Chose which arise in a 
variety of schemes. Of course, not every class 
would be natural Co every theory; and not every 
possible scheduling class can be encompassed by 
any one viewpoint. The idea, nonetheless, is to 
generate, as it were, some of the characteristic 
features encountered in scheduling. 

For the above mentioned, illustrative, pur- 
poses, then, we select Che Total Synthesis Problem 
(TSP) structure proposed by Peczkowskl, Sain, and 
Leake [5] in 1979. In the TSP idea, both the com- 
mand/output-response, represented say by a matrix 
[T(s)I, and the comniand/control-response, repre- 
sented say by a matrix [M(s)J, are to be simul- 
taneously synthesized subject to the constraint 
imposed by a plant, which could be represented by 
a matrix (P(s)]. Fundamental to TSP is the Nomi- 
nal Design Equation (NDE) 

(T(s)] - [P(s)HM(s)] , 

which must be satisfied no matter what controller 
structure might be selected. For a full discus- 
sion of NDE, see [6]. Once NDE is solved, one can 
make use of any of the multitudinous structural 
synthesis methods of modern theory or modern prac- 
tice CO develop a feedback controller. Feedback 
synthesis generates a number of interesting ques- 
tions of theory. For examples of these and for 
additional references, see [7]. It also generates 
very Important questions of practice. Examples of 
these have been given in (8,9,10). In [SJ , the 
structure of Figure 1 was employed on a turbojet 



Figure 1. Sample Controller Structure. 
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example. Subsequently, (9] and [10] considered 
particularizations of [8], especially insofar as 
L and K are concerned. 

The structure [8] of Figure 1, however, is 
adequate for purposes of subsequent discussion. 

Let farailies {L,K,G,H} have been determined by a 
technique of the reader's choice, in correspondence 
to a family {P} developed along desired operating 
lines of the plant. From a more global point of 
view, the plant may be regarded as a general, not 
necessarily linear, function p : U Y for appro- 
priate control and response spaces U and Y. We 
are Interested in considering four cases of sched- 
uling: 


(I) 

L 

art. 

K 

as 

a 

function 

of 

r; 

(II) 

L 

and 

K 

as 

a 

function 

of 

their output 

(III) 

G 

and 

H 

as 

a 

function 

of 

Lr; 

(IV) 

G 

and 

H 

as 

a 

function 

of 

y> 


Notice that, in I and II, we may as well choose 
either L or K, because the issues are the same. 

Abstract Derivatives 

As Indicated in the Simple Example Section 
above, the idea of polynomlc scheduling, of gains 
or time constants, suggests a state description in 
terms of series. Because we wish to use operator 
theoretic methods to some extent, it is convenient 
here to make a few introductory remarks about der- 
ivatives in such a context. 

Let V and W be normed real vector spaces, 
with Z open in V. A function f : Z -*■ W is 
differentiable at a point p in Z if there ex- 
ists a continuous linear map F : V -»■ W such that, 
for (p+h) in Z and h in V, 


If F exists, then it is unique and is called the 
derivative of f at p, and is denoted by 

(Df) (p) : V -*■ W . 

In case f is differentiable on Z, then we have 
a construction 

Df : Z -*• L(V,W) , 

where L(V,W) denotes the real vector space of R- 
linear maps 7 - W. Higher order derivatives are 
defined in a recursive fashion, 

(D"f)(p) = (D(D-'^f))(p) , 

with r a positive integer, provided that the in- 
dicated limit exists. For more discussion of 
these notions, the reader may wish to consult [11, 
12,13). 

An important connection exists between the 
calculus on normed vector spaces and the tensor 
algebra. Indeed, 



D^f(p) e L(V,L(V,W)) , 

D^f(p) e L(V.L(V,L(V,W))) 

whenever the limits exist. Let us denote by 
L(V^,V2,...,V^,W) 

the real vector space of n— linear functions 

V, K V- X ... X V W , 

1 Z n 

an n-llnear function being one which is linear in 
its remaining argument whenever (n-1) of its 
arguments are fixed. It can be shown that there 
exist Isomorphisms 

L(V^,V2,W) - L(Vj^,L(V2,W)) , 

L(V^,V2,V3,W) L(V^,L(V2,L(V3,W))) , 


so that (D^'f) (p) can be regarded as an r-linear 
map 7*^ ■+■ W, up to isomorphism. We suppress this 
Isomorphism and think of (D*^f)(p) as just such a 
map . 

It is now straightforward to establish a con- 
nection with the tensor algebra, and we do so in 
Che section following. The importance of the con- 
nection lias in its parametric possibilities: 

Every r-linear map can be composed from a linear 
map and a universal r-linear consti.ut.cion called 
Censor product. In a sense, the linear map em- 
bodies the parameters which are available for 
scheduling; and we pursue this view in a later 
section. 


Tensor Algebra 


In this section, we develop some of the struc- 
tures with which we can subsequently discuss sched- 
uling questions I-I7. Let 7 be a real vector 
space. For each integer r which is two or great- 
er, let 


(o’^7, a’^) 


be a Censor product for r copies of 7. The no- 
tion extends to 1 and 0 by Che definitions 

aH » 7 , a*^7 = R . 

T3ie sequence a’^7, r = 0,1,2,..., can be developed 
into a biproduct, and the Images of a^7 under in- 
sertion can be given the same notation. Then the 
tansorial powers can be developed into an as- 

sociative algebra by defining Che internal direct 
sura 

03 

= I a"v , 

n=0 
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and by 

equipping 

av 

with 

Che 

(0,0) -*• 00 for 

a, 0 

, O 0 

e 8 ' 

defined by 






00 - 

1 

0 8 




n,m 

(1 

where 

a ■ ) a , 
^ n 

6 ■ 

I 8 

t m 

for 


n ' m 

s a“v. with this multiplication, aV becoraes the 
graded tensor algebra over V with elements 


which are sequences of the tensors e 

a^V, 1 “ 0,1,..., and with unit element (1,0,...). 
We emphasize the fact that multiplication in the 
tensor algebra is not a tensor product. 

Now let ®V and aw be tensor algebras as 
defined above, over V and W respectively. For 
every pair n,m a 1, let a”(V,W) be a tensor 

product of a'V and a”V, that is. 


product ®2 on the right is between an (n+m)-ten- 

sor and an (r+a)-tensor. 

An algebra structure may be placed on ®(V,W) 
by defining a multiplication operation. To this 

end, let o” e ®”(V,W) and 6 e 8 (V,W) so that 
mm 3 s 

the tensors 


a 


I 

n^m 


n 


a 


m 


0 - I 

r,s 


are elements of ®(V,W). Then the product of two 
such censors is given by 

ID S 

n,m 

r,s 


where the symbol ® is the same as ^2 above. 
Notice that the multiplication rule Implies 


®”(V,W) - (®"v) 8 (a“w) . 


(o ® 0 )(a 8 0) 

n m r s 


(a ® 


®9 

m z r 


3s> 


We sec 8”(V,W) - ®*V and ®°(V,W) » a “ w . In a 
U m 

manner similar Co that preceding, 

a"(V,W) , n - 0,1,2 m - 0,1,2,... 

□1 

can also be developed into a blproducc; and Che 
images of each of these spaces under natural inser- 
tion into Che hlproduct can again be given the same 
symbolic representation. Again, then, we construct 
Che internal direct sum 

00 

0(V,W) = I ®"(V,W) 
n,mS0 

with 


« (o a o ) 8, (0 ® 0 ) 

n rim s 

- (a a ) 8 (0_0_) • 
nr ms 


This relation shows that the algebra ®(V,W) is 
the canonical tensor product of the subalgebras 
®V and SW, or 

8(V,W) - (3V) ® (8W) . 

The results in this section are adjustments 
of Chose which may be found in [14]. Our motiva- 
tion is, of course, the expansion of functions 
f : X X u X, for X a real vector space of 
states and U a real vector space of controls. 


®(V,W) » I [ I «"(V,W)] 
k “0 n+m»k 

functioning as the Induced gradation on ®( V , W ). 

Now consider two spaces ® (V,W) and ® (V,W). 

m s 

There exists a unique bilinear mapping 

U : ®"( V , W ) X ®’^( v , W ) ®" t !( V , W ) 

m s mrs 


Formal State Descriptions 

Consider a nonlinear state descr;i.pclon of the 
general form 


X - f(x,u) 


for 


f : X X u X 


with action 

lj(a ®8,a ®0) = (a 8a)® (S 80) 

nmrs nr ms 


with X and _U_ real vector spaces, equipped with 
norm. Let (x,u) be a point in X x u, and sup- 
pose Chat 


where a s 


a e ®”^V, B e ®^, 0 £ 

r m s 


s 


W . The 


pair (a , (V,W) ,p) is a censor product, or 

®!1^’^(V,W) = ®"(V,W) a a*^(V,W) 
mTS m s 

and 


(a„ 8 ap ®^ (0^ 8 0^) - (a^ 8 s^) 8^ (a^ 8 0^) 


We have subscripted the product symbol a in this 
equation In order to emphasize Che fact Chat the 
defining product ®^ on the left is between an 
(n+r)-censor and an (nri-s )- censor , while the defined 


D^f : Z f L(X X X u,x) 

is available for r = 0,1,2,..., with Z open in 
X X u and (x,u) in Z. Then, formally, 

f(x + x,u + u)» \ ^ (D^f) (x,u) (x,u) 

k =0 

(k) 

where (x,u)' =■ ((x,u) , (:c,u) , . . . , (x,u)) , the 

right member having (x,u) k times. It should 
be recognized that this series could be replaced by 
a finite number of terms together with a remainder. 
However, Che above representation is adequate for 
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brief iLtuscraClve purposes. Space does not permit 
a discussion of whether, or how, the series accept- 
ably describes the function. Along the same lines, 
we pass over the related question of how it affects 
the vector field associated with the differencial 
equation, and therefore its solutions. Instead, we 
remind Che reader that (D^f)(x,u) is a k-llnear 
function on (X x U)*^ to Xj and this suggests 
Chat we can use tensor algebra to parameterize it. 
Indeed, denote by (x,u)^ the k-fold tensor pro- 
duct of (x,u) with itself. Then Che k-llnear 
function (D^f)(x,u) can be factored uniquely in 
the manner 


(x,u) 9 (x,u) - l^(ir^(x,u) 9 ir^(x,u)) 

+ 1 (it (x,u) 9 ir (x,u)) 

UA U A 

If we identify Images of Che injections with their 
domains, as for example 

1 (U 9 U) - U 9 U , 
uu 


Ljj(x.u) 


Chen we can write 


where 


(x,u) 9 (x,u) “X9x + x9u + u9x + u9u . 


(9^(X X U), 9^ 

is a tensor product for k copies of X x u, or 
what is sometimes called a kth Censorial power for 
X X U. In this case, the kth parameter map oper- 
ates in the manner 

Lj^(x,n) : 0^(X X U) H- X . 

We have, therefore, r’’’’at 

f(x + X, u + u) " I L^(x,u) 9 9'^(x,u)^'^^ 

k*"0 

" ^ ^ L^(x,n)(x,u)*^ . 

k”0 


According to the conventions of 9(X,U), however, 
discussed in the section preceding, we agree to 
write 


U 9 X 


T X 9 u 

ux,xu 


for an appropriate Isomorphism 

way, we can proceed to 

- - 2 - - 2 
L2(x,u)(x,u) - l2(x,u) X 


+ I.2(x,u) X 9 u 

+ L„(x,u)T X 9 u 

2 ux,xu 


In Chat 


+ L2<x,u) u , 


Next consider the rearrangement of a term of 

type 

tm am 

L^(x,u)(x,u) . 

Consider, for example, the case k - 2, namely 
2 

(x,u) =• (x,u) 9 (x,u) . 

Such a form does not relate directly to the struc- 
ture of the section preceding, which would Involve 
terms of type xJ 9 u®. However, there is a nat- 
ural way to convert to that structure. Define pro- 
jections 


which we re-noCate to (with factorials included) 

__2 -- --2 

L 2 q(x,u) X + L^^(x,u) X 9 u + Lq 2 (x,u) u . 

In this way, the formal expansion becomes 

“ ® j , 

f(x + X, u + u) “ I ^ L. ,(x,u) X 9 u^ 
i-0 j=0 

from which point we can examine the scheduling 
questions previously raised. 

Internal Scheduling 



Tf : 

X X 

u 


u : 

IT : X X 

U * X ; 


u 




X 


and 

inj actions 







i : 

U 9 

u 

-w 

S ; 

i : U 

9 X -<■ S ; 


uu 





ux 



i : 

X 9 

u 

-h 

S ; 

i : X 

9 X -<■ S ; 


XU 





XX 


for 









S = (U 9 

U) > 


(U 

9 X) 

X (X 9 U) 

X (X 9 X) 


In the section on Multivariable Motivation, we 
brought attention to four cases of scheduling. With 
Che aid of the ideas foregoing, we would now like 
CO comment on each of these. We shall see chat the 
scheduling idea of our Simple Example Section, 
while motivating in nature, is not rich enough to 
embody the complete idea in question. In particu- 
lar, we have primary interest in the scheduling of 
parameters as a function of (x,u) and not (x,u). 
Though such issues would have encumbered the Simple 
Example, we have established now enough background 
to make the consideration. 


Then we can write 


Case 1 


It is sufficient to consider K as a function 
of r. Notice chat r is effectively u for this 
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case. Accordingly, we have as a mapping 

: U * L(0j(X,U),X) , 

and we may write 

CO 00 

fCx ■■ X, u + u) - I I L. ,(u) x^ ® . 

i-0 j-0 J 


Suppose, for purposes of G, that H is an iden- 
tity and that plant output is its state. What we 
can do is to examine the series combination of F 
and G, with state 

Z - (Xg,Xp) £ Xg K Xp . 

If we take u as the input to G, then 


The constraint Imposed by Case I is thus a domain 
restriction on form that the action 

of depends only upon u and not upon x. If 

we make the formal expansion 


Lij (u) 


I 

k-0 


r "I' 

^ijk 


in a manner analogous to the steps taken above for 
f : X X u X, then 

CO 09 00 

f(x + X, u + u) - I I I L u (x^ ® u^) . 
i-0 j-0 k-0 ■' 

Assuming distribution, we find that the structure 
of Interest Is engendered by 

Lljk “ ® u-^) . 

This raises the forms 

u ® X ® u-^ , 

which can be referred, by isomorphism, to the pre- 
viously described tensor algebra by 

i „ -k - j 
X 0 u ® u . 


Case II 

We shall assume, for simplicity, that y is 
equal to x. Again, it is sufficient to consider 
Just one, say K, of the two mappings. In a man- 
ner similar to Case I, it can be shown that 

f(X + X, U + U) “ I I [ ^<41, ® u^) , 

1=0 j=0 k=0 

witih the corresponding underlying construction 

-k 1 1 

X ® X © u*^ . 


z ■ fUfU) 


with f : Xg X Xp X U X Xp and with 

u • r *• TT 2 


for IT 


: Xq X Xp Xp an appropriate projection. 
Again, we have the basic expansion 


CO 00 


f(z + z, u + u) “ I )] L .(z,u) z^ ® u^ , 
i-0 j-0 ^ 


and we wish to schedule with 

r - u + IT z . 

We can reduce this situation to one which is close 
to chat of Case II. Write 


L., (z,u) - L. . (z, r - TT z) 


“IJ 


ij 


I I • 


k-0 m-0 


ijkm' 


Then the basic issues involve 

-m -k i j 
r © z © z © u , 

The remaining part of this case has to do with H. 
A similar approach can be applied, if we assume G 
to be an Identity. Here 

z “ ’ 

and 


u » r - TT z i 


Remark 

Although it has not been stated as one of our 
cases , Che generalization of Cases I and II clearly 
leads to 

-k k -k k 
X ® X ® u ® u 

Case III 

This case is quite a bit more complicated Chan 
Che previous two cases. For simplicity, cake L 
CO be Che identity and K to be zero. In consid- 
eration of G, we note chat its input, say u, Is 
a sum of r with a function of the plant output. 


the results vary essentially only in the meaning of 

z . 

Remark 

Assumptions on L, K, G, and H are for con- 
venience only. They can be removed easily by ex- 
panding z and the definition of input. 

Case IV 

Notice chat the question of H schedTiled on 
y is Che same as H scheduled on its input; and 
this reduces to Case I. Moreover, the question of 
G scheduled on y is essentially that of G 
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scheduled on a projection of Che sum state 

Xp) , which In cum is a question of the composite 

loop system scheduled on its state, which in cum 
reduces to a version of Case II. 

These are a few of the basic concepts which 
arise in nonlinear scheduling of state equations. 

We have introduced a viewpoint which seems broad 
enough to permit discussion. We have not used this 
framework to solve any new problems, however, but 
only examined its possibilities for providing 
frames of reference. 

Once internal scheduling has been carried out 
in an acceptable way, one comes next to Che re- 
lated inpuc/oucpuc functions. This is also a large 
subject, and we give only a sample survey of the is- 
sues in the next section. 

The Inpuc/Output Connection 

In this section we briefly highlight some gen- 
eral results from the literature relating Che in- 
ternal and external constructions for nonlinear dy- 
namical systems. Such results suggest just a few 
of Che approaches to the mathematical bridge nece- 
ssary for the scheduling discussions of the present 
paper. 

We consider a nonlinear control system 
X “ f(x,u) 
y ■ h(x) , 

where f :XxU-*-X, h!X-*-Y, for X, U and Y 
Che state, input and output spaces, respectively. 
With appropriately defined conditions, there is 
associated with such a system an input-output map 
from U to Y. Moreover, if we assume a suffi- 
cient degree of smoothness for f and h, this 
input-output map may taite on a functional expansion 
representation commonly known as a Volterra series, 
first considered by Vol;erra [15] as early as 1890, 
and studied later in d« pth by Weiner. 

An important subclass of such systems is that 
of bilinear systems, which in recent years has it- 
self motivated the study of Volterra series. In 
[16] Che construction of, and realization from, 
Volterra series representations for bilinear sys- 
tems is outlined, where it is shown that a neces- 
sary and sufficient condition for the existence of 
bilinear realizations corresponds to that of fac- 
torizabllity of the kernels which characterize the 
Volterra series expansion. In light of Krener's 
result [17] for the bilinearization of a rather 
general class of nonlinear systems, such realiza- 
tion studies '■ake on considerable utility. This 
latter result is based on a linearization technique 
used by Brockect [18], introduced first by Carleman 
[19] in 1932, and used later for similar applications 
by Bellman and Richardson [20] in 1963. We note 
here that it is a straightforward exercise to de- 
velop this technique analogously in terms of the 
tensor algebra and corresponding products. 

Brockett [21] has expanded these ideas to a 
more general class of systems for which 

X - g, (x) + [ u.g.,(x) 
i 

y “ h(x) , 
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where g^^, ^2 ' "*■ X and h : X Y are analy- 

tic, and Che u^ are components of the vector u 

€ U. Such systems are termed linear analytic and 
are closed under composition and feedback. It is 
shown in [21] r.nat the solutions can be expanded in 
a Volterra sex.ies provided that there is no finite 
escape time. Furthermore, necessary and sufficient 
conditions for n Volterra series to be realizable 
as a linear analytic system are given. For a fin- 
ite Volterra series representation, this will be 
the case if and only if the kernels are separable. 

In [22], Gilbert develops similar results for 
the general system and for the linear analytic sys- 
tem. The approach is an alternative to the Carleman 
technique, and uses Frcchec power series in the 
functional expansion. Convergence results of 
Brockect [21] are utilized for truncated Volterra 
series. 

Yet another approach to these results is in- 
troduced in [23] where standard tools from calculus 
and sralysis are employed in studying the existence 
and uniqueness of Volterra series representations 
for nonlinear systems. In their existence proof 
for the Volterra expansion of the class of linear 
analytic systems, Leslak and Krener exhibit a pro- 
cedure for constructing Che kernels of the expan- 
sion. This construction may be applied to more 
general systems, where the input enters nonllnear- 
ly, by first approximating the system by a bilin- 
ear system, using the result in [17]. In the work 
of Crouch [24], the rich mathematical structures 
of Lie groups and differencial manifolds are ex- 
plored in relation to the realizations of finite 
Volterra series for linear analytic systems, A co- 
ordinate free development is presented as an exten- 
sion to the results of [21] and [23]. 

Truucation 


The preceding remarks have made use of series 
representations for internal mcidels and have dis- 
cussed relations between Internal and Input/output 
models. It is, however, unllkexv Chat one would 
schedule an entire series in practice, unless chat 
series were simply an alternate form for a function 
which could be parameterized in finite terms. Thus, 
although the series format permits a rapid intro- 
duction to certain of the underlying structures in- 
volved in scheduling, we wish to mention some is- 
sues which arise when Che series are truncated, or 
when they had a finite number of terms from the 
outset. 

A principal feature of truncation in an in- 
ternal or state model is Che resulting approxima- 
tion which occurs in the input/output model. On 
appropriate normed real vector spaces U and Y 
of controls and responses, respectively, one might 
express Che input/output function by 

p : U -*■ Y 

and, formally at least, its series expansion by 

CO 

p(u + u) - ) d’‘p(u)u^''^ . 

k-0 

The algebra of input/output mappings, under compo- 
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sitlon and polntwiso addition, Is of ringlike na- 
ture, falling only one distributive law. Minor 
conpllcatlons arise, moreover, because U s* Y. 

It Is important to realize chat similar remarks 
can be made when the series representation is 
truncated, as for example In the manner 

p(u + u) - I — • 

k"0 

Call Che right member of this equation 
P^(u.u) . 

Then we have fundamental questions to be addressed 
Insofar as connections are concerned. 

Example 1 (Parallel) 

Let p^ : IJ -*• y, 1 » 1,2. Then under what 
conditions does 


Example 2 (Series) 

1 2 

Let p : U •* V and p : Y Z. Then un- 
der what conditions does 

,2 1 , 2 1 . 

(p op) “P °P ’ 

m “^m *^m 

Of course, the feedback case can also be addressed. 

Such questions have been addressed In [25,26,27]. 
Though we have no space remaining here for further 
detail, an important point Is that the questions 
can be addressed within the framework brought for- 
ward in this paper. This offers possibility, then, 
for determining features of scheduling which in- 
volve interconnections in the controller of G, H, 
L, and K. 

Conclusions 

Scheduling of local, linear multivariable con- 
trollers into global, nonlinear multivariable con- 
trollers is a very present tool In a variety of 
modem control applications. A case in point is 
gas turbine engine control. See, for example, the 
companion paper [4]. 

Despite the entre.'.chment of scheduling in 
practice, the theoretical view of such approaches 
is relatively undeveloped. There are many reasons 
for this. Including the fact that theory is pre- 
disposed to give general conditions across many 
classes of possible application and finds it qul.'-a 
challenging, for the most part, to be general about 
particular applications. Nonetheless, it seems 
that much could be gained by a theory/design dialog 
on the issue of scheduling. 

Toward such dialog, we propose in this paper 
a viewpoint which may be broad enough and deep 
enough to support such a dialog. The central fea- 
ture in the viewpoint is the parameterization pos- 
sibility which is inherent in tensor algebra. For 
a reasonably general control configuration, we 
have used the tensor framework to examine four 
typical scheduling structures. Though much remains 


to be done in such a dialog, we have found these 
ideas useful in thinking about the Joint theory/ 
application issues which are involved. 
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A COMPUTER-AIDED DESIGN PACKAGE 
FOR NONLINEAR MODEL APPLICATIONS* 

T.A, Klingler, S. Vurkovich, and M.K. Sain 

Deparcmenc of Electrical Engineering, University of Notre Dame, 
Notre Dome, Indiana 46236 


Abstract . An important aspect of multivariable control system design in- 
volves the formulation of reliable mathematical models. Gas turbine engine 
control systems, with their inherent nonlinearities, provide common practic- 
al examples of the need for nonlinear models. In this paper we discuss a 
computer-aided design package for generation of such nonlinear models, using 
an approach involving notions of power series and algeoraic tensors. Two 
independent computing systems are employed Interactively in the overall pro- 
cess of model formulation, identification, and validation. The package is 
sufficiently generalized to accommodate any particular nonlinear modeling 
problem when formulated within the framework of the algebraic tensor scheme. 

Keywords . Computer-aided system design; multivariable control systems; mod- 
elling: tensor algebra; nonlinear systems; algebraic system theory. 


INTRODUCTION 

.Models have always been an Important aspect 
of applications engineering in the area of 
multivariable control system design. See 
for example the work of Kreindler and 
Rothchild (1976). Practical and industrial 
examples of the use of models are provided 
by gas turbine engine control systems, which 
commonly use models to generate control and 
response trajectories for various power de- 
mands. These models, when scheduled over 
the operating envelope, can reduce the com- 
pensation normally required of the control- 
ler, and thus provide the feedback loop with 
an opportunity to achieve better accuracy in 
the presence of noise and parametric uncer- 
tainties . 

The scheduling of feedforward models and 
feedback compensation typically produces 
nonlinearities, even if the local models are 
linear. Accordingly, there is basic inter- 
est in fundamental approaches which incor- 
porate nonlinearity at the outset. Such ap- 
proaches should (1) reduce to the earlier 
linear schemes for variables with small ex- 
cursions, (2) be amenable to scheduling, and 
(3) offer opportunities for determination of 
parameters from simulation data. 

One such approach, investigated by Yurkovlch 
and Sain (1980) and Klingler, Yurkovlch, and 
Sain (1902) , uses the notions of power ser- 
ies and algebraic tensors (Sain, 1976) to 
generate a class of nonlinear models. The 
important feature of the algebraic censor is 
that it provides an organized way of de- 
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scribing the power series expansion formula, 
lending itself with relative ease to pro- 
grannnlng on a digital computer. Furthermore, 
its use allows for the implementation of 
linear parameter identification techniques. 

This paper reports on Che develonment of an 
interactive computer-aided design package 
for the formulation, identification, and 
validation of one particular model structure 
which uses Che above-mentioned tensor appro- 
ach. The software package, developed using 
Che extensive capabilities of the IBM and 
DEC Command Procedure Languages with the 
strengths of FORTRAN and SPEAKEASY, is di- 
vided into two segments and tailored to 
utilize existing computer hardware affec- 
tively, as well as to provide the fastest 
possible user turnaround time. The inter- 
active nonlinear model generation segment 
is lEwlemented on a Time Sharing Option 
(TSO) of Che IBM 370-168 computer system, 
where the memory dependent and highly com- 
putational routines of the package can bene- 
fit from use of the virtual memory and float- 
ing point hardware. Once a structured non- 
linear model is identified, it is then 
tranaferev' to the DEC PDF 11/60, where in an 
interactive environment it can be analyzi’d 
and compared to a linear model as well as 
the true system. In this manner, the user 
has at his disposal both the graphics capa- 
bilities of the video terminal and an elec- 
trostatic prlnter/plotcer for the immediate 
display of data and comparative trajectories. 

The remainder of the paper is outlined as 
follows. First, we briefly discuss notions 
from analysis and algebra which form the 
foundation for the tensor approach used In 
the model formulation. A detailed discus- 
sion of the interactive design package is 
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Chen given, followed by a brief discussion 
of Che coiapucacional aspects regarding 
floating point operations in the model sim- 
ulation phase. We close with an e.xample 
problem from a turbofar. Jet engine simula- 
tion. 


XONLINEAR MODEL FORMULATIOH 

Prior to proceeding to the description of 
the computer-aided design procedure in the 
modeling scheme, we outline here some of 
the prerequisite mathematical Issues in a 
coordinate-free development. Since the 
treatment is brief, the reader may wish to 
consult Dteudonne (I960) and Greub (1967) 
for complete expositions of the topics dis- 
cussed herein. 

Tensor Ideas 


We begin with a discussion of abstract der- 
ivatives and the calculus of normed vector 
spaces. Let V and W be normed vector 
spaces and let Z be open in V. Suppose 
chat f : Z W is differentiable at a 
fixed point p in Z. Then the derivative 
of f : Z - W at p is a mapping 

(Df)(p) : V - W (1) 

where 

Df : Z -<• L(V,W) ; (2) 

chat is, the derivative mapping in Eq. (1) 
is an element of the real vector space of 
linear mappings from V to W. Higher or- 
der derivatives are defined recursively as 

(D‘^f)(p) - (D(D'^"^f))(p) . (3) 

for the positive integer r, provided that 
the (r-l)st derivative is differentiable. 
Moreover, higher order derivatives are them- 
selves linear mappings according to 

D^f(p) e L(V,L(V,W)) , 

D^f(p) £ L(V,L(V,L(V,W))) . (4) 


If L(v'^;W) denotes the real vector space 
of r-linear mappings from to W, it 

can be shown Chat there exist isomorphisms 

L(V^|W) -*• L(V,L(V,W)) , 

L(V^W) - L(V,L(V,L(V,W))) , (5) 


so Chat Che rch derivative of f at p can 
be regarded as a mapping from to W. 

We suppress this isomorphism and consider 
D"f(p) as an element of L(V'^;W). 

We now use this mulcilinearicy of the deriv- 
ative mapping to make a connection with no- 
tions from algebraic tensors. Let 
be a tensor product for r copies of V. 
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Recall chat by the unique factorization pro- 
perty of the censor product, for every map- 
ping i() j * W in L(V*^;W) there exists 
a mapping \ : ®*’V • W in L(»‘^V,W) such 
chat ij; <■ 1 • 9*^ for 9^^ : •* 9 V in 

L(V*';9*'V), Furchormore, the implication of 
Che unique factorization property is chat 

L(®’^V,W) L(V*^;W) (6) 

is a vector space isomorphism. Thus, via 
the Isomorphisms of Eqs. (3) and (6), Che 
r-linear mapping D^’f^p) : -*• W can be 

composed from a linear mapping W and 

the universal r-linear censor product map- 
ping 9'=’ : - a^v. 

State Description 


The ideas discussed above are now used Co 
formulate Che model structure for a given 
nonlinear system. We consider systems whose 
states and inputs are elements of the normed 
real vector spaces X and U, respectively, 
and which may be described by the nonlinear 
ordinary differencial equation 

k - f(x,u) (7) 

for f : X X U - X. Let (x,u) be a fixed 
point in Z open in X x u, and suppose 
that f : X X u ♦ X is of sufficient smooth- 
ness on Z. Formally, 

ao 

f(x+x,u+u) " [ rf (D^f ) (x,u) (x ■ i) , 

k"0 ’ i ) 

(k) 

where (x,u)^ - ( (x,u) , (x,u) , . . . , (x,u) ; k 

times. Due to space limitations we cannot 
address existence or convergence questions 
relative to Eq. (8). We note, however, that 
this series could be replaced by a finite 
number of terms together with a remainder 
term in a standard application of Taylor's 
formula. 

•According to chn preceding discussions we 
now make Che unlviue factorization 

o'^f(x,u) - Lj^(x,u) <■ , (9) 

where L. (x,u) : ®^(X x tJ) - x is linear. 
Denote the k-fold censor product of (x,u) 
with itself by (x,u)‘^ so chat, with Eq . 

(9) , we have 

00 

f(J^x,u+u) » I ~ Lj^(.x,u) (x,u)''. (10) 
k”0 

Sain and Yurkovich (1982) have shown that 
Che individual terms in the series of Eq . 

(10) may be rewritten as, for e.xample in the 
case of k " 2, 

1 - - 2 

jf- LjCx.uKx.u) - Ljq X 9 X + 

L^j^ X 9 u + Lq 2 u S u , (11) 

where we have suppressed the notation of 
(x,u) on the right side of Eq. (•’ ■'. In 
this way Eq, (8) becomes 


(M 01 

f(x+x,u+u) - J* V L. a u , 

i-O J-O J (12) 

forming tho structure for the nonlinear mod- 
el. 

Application 

In practical appilcationa a truncation ap- 
proximation of tho series in Eq. (12) la 
conaidared. In corma of computing, then, 
the teak in Che modal building schomo 1« to 
identify the paramocora concslned in matrix 
reprauontatlona of the Lj^,. Ordered buses 
in K, U, X 9 ,K, .K s) U, and ao on, are 
cl\oson a priori to be used ■ .'uacructing 
a linear louse squares i Im ‘ . ,, cion prob- 

lem. Tho i^rdorlng aigori' ■ i v.arkovlch, 
1981) which faclllcatea this procedure, 
used in tho incoracclvo software package de- 
scribed heroin, is easily impiomoncod on u 
digital compuCor. 

In practice, a differencial equation doscri- 
pcion of tho nonlinear system may or may not 
bo available. In either case, the basic 
strategy involves initial condition and con- 
trol signal design so that, through data 
sampling and derivative oacimution, a model 
of the original svstom of Eq. I?) may bo 
identified. Tho nonlinear model la required 
not only to outperfonn a stand.srd linear 
model locally about an expansion point, but 
to establish a larger region of model val- 
idity as well. 

CATN/VP 

Tho Intent of this section is to present a 
detailed discussion of the Computor-Alded 
_renaor 2<onlinoar Modeling Applications pack- 
age iCATNAP) currontiv used as a dovolopmant 
tool in tho formulation, Idontlf ication, and 
validation of nonlinear models of the type 
mentioned above. The structure of CATNAP 
is bused upon Ideas from distributed proces- 
sing and Local networking (Tanenbaum, 1981) 
In which computations are spread over mul- 
tiple machines. More specif lc,sllv CATN,U’ 
is divided into two segments, each of which 
is implemented on an IndeponTlont computing 
system. These segments are entitled CENER- 
,ATE and SIMUL.\TE. CENr.R,\TE is Inip lamented 
on a Time Sharing Oicion ^IBM, 1975) of tho 
IBM 370-1 p 8 mainframe computer and is used 
to formulate and identify models, whereas 
SIMUIATE is Implomontad on tho DEC PDP 11/ hO 
computer system and is used to study model 
validity and performance. Furthermore, both 
of these segments are highly interactive and 
contain straightforward input prompts as 
well as informative error mossagos, 

CL.,~ERi\TE 


Tho CENER/VTE segment of CATNAP Ij prlmarllv 
made up of three routines governed bv a 
higher level supervisor. Figure 1 contains 
a block dliigram depicting the structure of 
CENERiVTE. 


original page 
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CENER^VTE Suporviaoiy . This suporvisor>’ 
level is written using the command procedure 
language CLIST (IBM, 1976) and performs the 
following main functions in sequence t 


(1) prompts the user for tho name 
of the desired loader routine 
to b'' executed; 

(2) passes control to tho loader 
routine doflnori in (1); 

(3) passes control to IDENTIFY; 
and , 

l4) upon user raquest, passes con- 
trol to TRilNSFER. 


In ad.iiiion to these main functions, cer- 
tain laainconance roles such as file crea- 
tion, allocution and deiorlon are handled 
by this supervisor. 



Fig. 1. Block diagram for the CATNAP 
segment GENERATE. 


Loader Routine . Asaoclatad with each non- 
linear svstom to be modeled, chore exlaca a 
loader routine which performs cho modal 
formulation cask, The.se routines are 
scored in a library and are typically writ- 
ten in double precision FORTRAN. 


Tho purpose of any loader routine la to e.x- 
clte the given nonlinear svstem via Initial 
condition and control Input perturbations 
and to sample cho acatos, inputs and deriv- 
ative ascimacoa over h selected points in 
time. The svstom is then represented bv the 
m.icrlx equation 


nxh 


f'-ui hn ^20 ‘'u ‘■30 

nxp 


! X„ 
pxh 
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The first n+m rows of Che matrix are 
formed from the sampled values of x and 
ui Che remaining p-(n+m) rows are formed 
according to the ordering algorithm prev- 
iously mentioned, which minimizes the num- 
ber of computations. X i-j formed by load- 
ing derivative estimates for xi,X 2 , • ■ . ,Xjj 
at the h time points. The number p de- 
pends on n and m, and the degree of the 
truncation approximation. All this data is 
then scored in TEMPFILE for later use. 

Using this approach, CATNAP can accommodate 
any particular nonlinear modeling problem 
since the problem specifics are transparent 
to the remainder of che package. The only 
requirement is that TEMPFILE contains the 
appropriate data. 

IDENTIFY . After the completion of any 
chosen loader routine, the program IDENTIFY 
is executed. IDENTIFY reads the interim 
data from TEMPFILE and forms a least squares 
minimization problem which is solved for the 
partitioned matrix containing the desired 

parameters. These L^i^j parameters are 
recorded at the terminal as well as entered 
into the MODEL data file. 

It should be noted here that IDENTIFY is 
written in the high level language SPEAKEASY, 
which is based on the concepts of arrays and 
matrices and processes these as entities. 

This results in the elimination of the many 
loops necessary in other programming lan- 
guages. See the work of Cohen and Pleper 
(1979) , The main reason for employing 
SPE.AKEASY here is chat the highly efficient 
routine SBIEQUAT can be easily used to solve 
the least squares problem via singular value 
decomposition, thus reducing the apparent 
complexity of the problem to a minimum, 

TRANSFER . Upon a yes response to a super- 
visory prompt, Che program TRANSFER is sub- 
mitted batch to Che IBM 370-168. TRANSFER 
is merely a Job Control Language (JCL) deck 
which sends a copy of the file MODEL, con- 
taining the parameters, to the DEC 

PDP 11/60 computing system by the way of a 
Remote Job Entry port, and stores it in Che 
nonlinear model library. An excellent ac- 
count of JCL can be found in Brown (1977) . 

SIMULATE 


(3) execute VERSATEC which pro- 
duces hardcopy plots. 

As earlier, thi.s supervisor performs a num- 
ber of file maintenance duties in addition 
to the above functions . 



Fig. 2. Block diagram for Che CATNAP 
segment SIMULATE. 


Simulator Routine . A FORTRAN simulator 
routine usually exists for each nonlinear 
modeling problem studied; however, only one 
subroutine in that program is altered among 
versions, and that is the application sub- 
routine TRUES. The remainder of the pro- 
gram stays unchanged. TRUES contains the 
true system representation of the nonlinear 
system being modeled, and is used exten- 
sively in comparison studies. Because of 
Che number of TRUES subroutines that exist, 
a library has been created to store the 
various simulator routines. 

The execution of a particular simulator 
routine can be divided into three steps; 

(1) problem configuration; (2) systems in- 
tegration; and (3) solution display. 


Shifting our concern away from the discus- 
sion of GENERATE, we now focus our atten- 
tion on the SIMULATE segment of CATNAP. 
Basically, two routines plus a supervisor 
comprise the structure of SIMULATE. Figure 
2 offers an illustration of this structure 
to supplement the following presentation. 

SIMULATE Supervisory . Written in Che form 
of an Indirect Command File (DEC, 1979), 
this supervisor allows the user to; 


The first of these steps requires the user 
CO decide which of the available systems, 
true solution, linear model and/or nonlinear 
model, should be included in the session 
configuration. When a model is chosen, the 
user is asked to enter the name of the de- 
sired model. That model is then read into 
Che simulator from Che appropriate library. 
The linear models used in CATNAP are gen- 
erally identified by standard techniques 
and are available for use in comparison 
studies . 


(1) create new simulator routines; 

(2) execute existing simulator 
routines; and. 


Next, the user is prompted for various in- 
tegration parameters such as scepsize and 
upper time limit as well as initial condi- 
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cions. Input amplitudes and frequencies. 

The configured systems are then integrated 
and the solutions are sampled at 100 points, 
evenly spaced in time. 

Finally, to assist in the data analysis, a 
number of options are available to the user. 
They include: 

(1) printing the solutions on 
Che Versatec; 

(2) displaying Che comparative 
trajectories on the Tektronix 
graphics terminal: 

(3) writing the trajectory solu- 
tions CO SPOOL for hardcopy 
plotting at a later time; and, 

(4) solving the configuration for 
another set of initial condi- 
tions and control inputs. 


The use of these options provides a powerful 
yet flexible capability for the study of 
model performance and validity. Further- 
more, when all three systems are included 
in Che configuration, an additional error 
criterion is generated and used in testing 
the performance of the nonlinear model ver- 
sus that of the linear model. 


V 

Let denote the absolute maximum error 

in the nonlinear model solution, as compared 
CO Che true simulation solution, over the 
time range of simulation for the ich state 
variable. Similarly, we define sj; for the 
linear model error. Then is the com- 
parison sf - sj;. Tims, if is negative, 

Che nonlinear model has exhibited a smaller 
maximum absolute error in the ith state, and 
in chat sense has outperformed the linear 
model. 


VERSATEC . The routine VERSATEC, written in 
FORTRAN, reads Che trajectory solutions from 
SPOOL and records at the Versatec printer/ 
plotter, a data sheet corresponding to each 
plot sec which follows. The comparative 
trajectories themselves are then plotted. 

MODEL SIMULATION 


:x • Lz , 


where x is Che n-veccor of states, L the 
parameter matrix, and z the p-vector of 
ordered monomial terms derived from the 
various symmetric products of x and u, 
the m-vector of inputs (Yurkovich and Sain, 
1980) . The least number of multiplications 
required to construct z is merely p-n+m, 
or the total number of terms which involve 
products. This number is given by 

p-(n+m) - 


i-1 

+ c 1 


j»l 


,n+(i-j)-l. -ari-J 
^ (i-j) j 



where d is the model degree. 


Assuming the use of a fourth-order integra- 
tion routine, Che number of FLOPs necessary 
to formulate and integrate the system as em- 
bodied by Che model, across one integration 
time step, is 4(n)(p). As an illustration 
consider a four-state, three- input model. ^ 
Suppose, for simplicity. Chat 100 integra- 
tion time steps is the equivalent of one 
second in real time. This translates 
roughly to 0.25 million FLOPs per second 
for a degree-3 model (an approximation which 
retains terms up to and including the third 
degree) . While there are many ocher obvious 
considerations involved in real time simu- 
lation, this number is well within the 
bounds dictated by state-of-the-art compu- 
tation speeds of ten million FLOPs per 
second. 


EXAMPLE 


In the example to follow attention will cen- 
ter around NASA's QCSEE ("Quixie") £uiec, 

£lean, £horthaul Experimental Engine. Wise 
(1974) provides an excellent overview of Che 
QCSE engine program. QCSEE is designed 
specifically for powerec-lifc , short-haul 
aircraft, and incorporates several new con- 
cepts not all currently used on turbofans 
to achieve operational efficiency in a 
quiet, clean manner. 


In this section we comment on the efficiency 
of the model structure discussed above by 
studying the number of floating point oper- 
ations (flops) necessary in typical simu- 
lations. It is common practice in computer 
architecture to design processors which re- 
quire no extra time for floating point ad- 
ditions calculated simultaneously with mul- 
tiplications. Thus, we concern ourselves 
primarily with the latter, and by FLOPs we 
will mean multiplies. Since the largest 
burden of Che computer In the simulation 
process is Che actual numerical integration 
of model differencial equations, we will 
analyze only that portion of the simulation. 

The system to be considered cakes the form 
of Eq. (13), or 


QCSEE APPLICATION 

For this nonlinear modeling problem, a com- 
plex eight-state, three-control digital 
simulation of the QCSE engine is employed 
(Mihaloew, 1981) . Using this digital deck 
as Che system representation, a loader rou- 
tine, QCSELOAD, is constructed to formulate 
a reduced order four-state, three-control 
analytical model. The engine states chosen 
are the combustor discharge pressure (P4GS), 
the core nozzle pressure (P8GS) , the fan 
speed (NL) , and the compressor speed (NH) . 
The control inputs used are Che fuel flow 

^This represents a typical model as inves- 
tigated by Klingler, Yurkovich, and Sain 
(1982) . 
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(WFM) , Che exhauac nozzle area (A18) and che 
fan pitch angle (BETAF) . In a similar way, 
che slmulacor roucine QIXSIM is builc using 
QCSEE as Che true system in che subroutine 
TRUES . 

For Che results presented herein, two models 
have been formulated using QCSELOAD at 922 
power demand: a second-degree nonlinear 

modal, and a first-degree linear model. 

Both formulations are made using 12 steady 
state perturbations in che state and control 
variables. Furthermore, the control Inputs 
are manipulated so that cosinusoidal signals 
are inserted. The observed states and in- 
puts are sampled over a six second interval, 
and the difference between these values and 
the corresponding equilibrium values, to- 
gether with che ordered censor product terms 
and state derivative values comprise che 
data necessary for the identification. The 
model parameters are easily computed from 
IDENTIFY and Chen sent to Che PDF 11/60 via 
TRANSFER. 

Using Che capabilities of QIXSIM and VERSA- 
TEC, several validation studies have been 
completed to dace, all yielding satisfactory 
results. Figures 3-7 contain a representa- 
tive plot sec from VERSATEC illustrating Che 
model performance for a particular input 
sec, as well as Che graphical capabilities 
of CATNAP. Table 1 contains a variable led- 
ger for Figures 3-7. 

CONCLUSION 

The importance of nonlinear modeling in mul- 
tivariable control systems could hardly be 
overemphasized. And Che applications side 
of che problem has benefited greatly with 
che advent of expanded and more versatile 
computing environments . 

Rarely does it happen, chough, that one com- 
puting system can accommodate all require- 
ments placed on it, particularly when pla- 
gued by multiple users demanding unlimited 
access. It often haopens , however, chat Che 
capabilities of several computing systems 
are at ones disposal, each with various fea- 
tures to offer. In this case schemes em- 
ploying Che notions of distributed proces- 
sing and local networking become extremely 
useful. 

We have presented one such scheme in che 
form of an interactive computer-aided de- 
sign package for a specific nonlinear model- 
ing problem. The software package facili- 
tates Che analysis of complex problems, with 
relative ease to Che user, from the initial 
model formulation and identification stage 
through CO the model testing and validation 
studies. Series ideas and algebraic censors 
are che main vehicles in che model formula- 
tion. The importance of che censor approach 
lies in its parametric possibilities, and 
ongoing research is currently underway to 
exploit further che richness of such struc- 
tures . 
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TABLE 1 Variable Ledger for Figures 3-7 


o : True Engine Response 

4 : Linear Model Response 

* : Nonlinear 

Model Response 

X., ; P4GS (psi) 

u : WFM (lb„/hr) 

X,: P8GS (psi) 

u,: A18 (in*^) 

x“: NL (rpm) 

u,: BETAF (deg.) 

x^: NH (rpm) 
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PROBLEM SUMMARY •••»•••••••• 

COHF IGURATlONi TRUE, LINEAR i NONLINEAR 

NUMBER OF STATESi 4 

NUMBER OF CONTROLSi 3 

LENGTH OF TENSOR TERM VECTOR! 33 

DEGREE OF APPROXIMATION! 2 


SOLUTION PARAMETERS! 


INTEGRATION STEPSIZEi 0.010 

UPPER TIME LIMIT OF INTEGRATION! 2.000 

NUMBER OF PLOT POINTS! 100 

SPACING BETWEEN PLOT POINTS! 0.020 


STATE NUMBER 


INITIAL CONDITION 

• 

• 

• 

ERROR criterion 

1 


-0.075 

» 

• 

* 

-0.533E*00 

2 


0.001 

• 

-0.9ISE-02 

3 


-0.250 

« 

m 

-0.698E*0I 

4 


0.500 

• 

• 

-0. 1 I IE*02 


CONTROL NUMBER 

• 

* 

■ 

AMPLITUDE 

• 

• FREQUENCY (CYCLES/SEC) 

• 

1 

• 

• 

74.290 

• 

• 

2.100 

2 

• 

-21 . 200 

• 

! .000 

3 

« 

• 

-0,239 

« 

B 

1 .500 

Fig. 3. 

Sample 

data sheet for 

the qcSEE 

example 



Fig. 4. Comparative solutions: Fig. 3, state 1. 
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Fig. 5. Comparative solutions: Fig, 3, state 2. 



Fig. 6. Comparative solutions: Fig. 3, state 3. 
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SOFTWARE DESCRIPTION FOR SECTION IV 

In this appendix, we provide the necessary documentation of the 
software that was used for both the first and second order controller 
problems. The computer programs are written in the high-level language 
SPEAKEZY [23] for use on the IBM 370/168 at the University of Notre Dame. 
The first program is FIRORDA, which provides the first order analysis of 
the example problem. This procedure, as was demonstrated, is identical 
to the solution of the Riccati equation in matrix form and substitution 
of the optimal control into the system equations. The program first in- 
tegrates the Riccati equation using the modified Euler method described 
in [28]. This provides the solution for the optimal value term V^Ct) 
and subsequently the first term in the controller expansion K^(t). A 
fourth-order Runge-Kutta integration routine then solves for the state 
variable x(t) after appropriate substitutions have been made for the 
controller term. The final part of the program provides for the reduc- 
tion of the array size (if necessary) and the plotting of both the regu- 
lated state variables and also the control variables. 

The second order analysis of the example problem required two pro- 
grams. The first program, called COMATRCS, provided for the calculation 
of the coefficient matrices that appeared in equations (4.4) and (4.5). 

In particular, the program calculated those coefficient matrices that re- 
main constant over the integration interval. These matrices may be listed 
as follows: 
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CMl - (A^ ® ® 


CM2 " (Am ® 

2.. 2 
X ' ’21 


CM3 


'02 ~ x' 12 '”x 

(A^, a + (1„ a A, 


'01 


X 


02 ' 

01 


“ ^10 ^ ^10 ^ ""lO ’ 

CM5 - + (1^ 8 Ag^) 

■CM6 - (Ajo 9 1^) + (1^ ® Ajg) . 


S?^ 


1 2 
1 )sr, + (1 
x“^ 21 ' X 


Aoi> , 


The program that actually does the second order analysis is called SECORDA. 
This routine integrates to get V^Ct), again using the modified Euler 
method, and subsequently calculates K 2 (t). Using the results from the 
FIRORDA program, the second order controller may be generated. The sys- 
tem is then integrated to yield the state variables by using a fourth 


order Runge-Kutta integration routine. The results are then plotted in 

the final section of the program. In the solution of equation (4.4) for 

K 2 (t) in terms of 7^(1) , we define 

VMl = 1 ® TT^ (K ® 1 ) , 

X s 1 u 

which needs to be updated as Kj^(t) varies. Also needed are the Kronecker 
product matrices 


El 

= 

1 ^ 
X 

® 

E2 

s= 

1 

X 

8 

S • 

E3 

= 


8 

8 

E4 

= 

’'l 

8 

K 2 . 

E5 

= 

Kj 

8 

Ki . 


Like the VMl terra, these matrices also must be updated at each integra- 
tion step. 
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The final section of this appendix lists the time-varying control 
gains and the optimal value expression terras. In each case, the a.^t-es- 

sions are listed as functions of time, which appears in the first column. 

The array KA lists K^(t); the array K2 a lists K 2 (t); the array 

VA lists V 2 (t); and the array V3A lists V^Ct). In all cases, we 

have assumed an integration interval of 0 to 5 seconds and a stepsize of 


pAgb Is 


0.05 seconds. 
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